AdaptiveTimeStepping_impl.hpp
Go to the documentation of this file.
1/*
2 Copyright 2024 Equinor ASA.
3
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18*/
19
20#ifndef OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
21#define OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
22
23// Improve IDE experience
24#ifndef OPM_ADAPTIVE_TIME_STEPPING_HPP
25#include <config.h>
28#endif
29
30#include <dune/common/timer.hh>
31#include <dune/istl/istlexception.hh>
32
33#include <opm/common/Exceptions.hpp>
34#include <opm/common/ErrorMacros.hpp>
35#include <opm/common/OpmLog/OpmLog.hpp>
36#include <opm/common/TimingMacros.hpp>
37
38#include <opm/grid/utility/StopWatch.hpp>
39
40#include <opm/input/eclipse/Schedule/Tuning.hpp>
41
42#include <opm/input/eclipse/Units/Units.hpp>
43#include <opm/input/eclipse/Units/UnitSystem.hpp>
44
46
48
49#include <algorithm>
50#include <cassert>
51#include <cmath>
52#include <sstream>
53#include <stdexcept>
54
55#include <boost/date_time/posix_time/posix_time.hpp>
56#include <fmt/format.h>
57#include <fmt/ranges.h>
58namespace Opm {
59/*********************************************
60 * Public methods of AdaptiveTimeStepping
61 * ******************************************/
62
63
65template<class TypeTag>
67AdaptiveTimeStepping(const UnitSystem& unit_system,
68 const SimulatorReport& report,
69 const double max_next_tstep,
70 const bool terminal_output
71)
72 : time_step_control_{}
73 , restart_factor_{Parameters::Get<Parameters::SolverRestartFactor<Scalar>>()} // 0.33
74 , growth_factor_{Parameters::Get<Parameters::SolverGrowthFactor<Scalar>>()} // 2.0
75 , max_growth_{Parameters::Get<Parameters::SolverMaxGrowth<Scalar>>()} // 3.0
76 , max_time_step_{
77 Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60} // 365.25
78 , min_time_step_{
79 unit_system.to_si(UnitSystem::measure::time,
80 Parameters::Get<Parameters::SolverMinTimeStep<Scalar>>())} // 1e-12;
81 , ignore_convergence_failure_{
82 Parameters::Get<Parameters::SolverContinueOnConvergenceFailure>()} // false;
83 , solver_restart_max_{Parameters::Get<Parameters::SolverMaxRestarts>()} // 10
84 , solver_verbose_{Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminal_output} // 2
85 , timestep_verbose_{Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminal_output} // 2
86 , suggested_next_timestep_{
87 (max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>()
88 : max_next_tstep) * 24 * 60 * 60} // 1.0
89 , full_timestep_initially_{Parameters::Get<Parameters::FullTimeStepInitially>()} // false
90 , timestep_after_event_{
91 Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60} // 1e30
92 , use_newton_iteration_{false}
93 , min_time_step_before_shutting_problematic_wells_{
94 Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day}
95 , report_(report)
96{
97 init_(unit_system);
98}
99
106template<class TypeTag>
108AdaptiveTimeStepping(double max_next_tstep,
109 const Tuning& tuning,
110 const UnitSystem& unit_system,
111 const SimulatorReport& report,
112 const bool terminal_output
113)
114 : time_step_control_{}
115 , restart_factor_{tuning.TSFCNV}
116 , growth_factor_{tuning.TFDIFF}
117 , max_growth_{tuning.TSFMAX}
118 , max_time_step_{tuning.TSMAXZ} // 365.25
119 , min_time_step_{tuning.TSMINZ} // 0.1;
120 , ignore_convergence_failure_{true}
121 , solver_restart_max_{Parameters::Get<Parameters::SolverMaxRestarts>()} // 10
122 , solver_verbose_{Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminal_output} // 2
123 , timestep_verbose_{Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminal_output} // 2
124 , suggested_next_timestep_{
125 max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() * 24 * 60 * 60
126 : max_next_tstep} // 1.0
127 , full_timestep_initially_{Parameters::Get<Parameters::FullTimeStepInitially>()} // false
128 , timestep_after_event_{tuning.TMAXWC} // 1e30
129 , use_newton_iteration_{false}
130 , min_time_step_before_shutting_problematic_wells_{
131 Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day}
132 , report_(report)
133{
134 init_(unit_system);
135}
136
137template<class TypeTag>
138bool
141{
142 if (this->time_step_control_type_ != rhs.time_step_control_type_ ||
143 (this->time_step_control_ && !rhs.time_step_control_) ||
144 (!this->time_step_control_ && rhs.time_step_control_)) {
145 return false;
146 }
147
148 bool result = false;
149 switch (this->time_step_control_type_) {
151 result = castAndComp<HardcodedTimeStepControl>(rhs);
152 break;
154 result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
155 break;
157 result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
158 break;
160 result = castAndComp<PIDTimeStepControl>(rhs);
161 break;
163 result = castAndComp<General3rdOrderController>(rhs);
164 break;
165 }
166
167 return result &&
168 this->restart_factor_ == rhs.restart_factor_ &&
169 this->growth_factor_ == rhs.growth_factor_ &&
170 this->max_growth_ == rhs.max_growth_ &&
171 this->max_time_step_ == rhs.max_time_step_ &&
172 this->min_time_step_ == rhs.min_time_step_ &&
173 this->ignore_convergence_failure_ == rhs.ignore_convergence_failure_ &&
174 this->solver_restart_max_== rhs.solver_restart_max_ &&
175 this->solver_verbose_ == rhs.solver_verbose_ &&
176 this->full_timestep_initially_ == rhs.full_timestep_initially_ &&
177 this->timestep_after_event_ == rhs.timestep_after_event_ &&
178 this->use_newton_iteration_ == rhs.use_newton_iteration_ &&
179 this->min_time_step_before_shutting_problematic_wells_ ==
181}
182
183template<class TypeTag>
184void
187{
188 registerEclTimeSteppingParameters<Scalar>();
190}
191
192// See Doxygen comment on the declaration in AdaptiveTimeStepping.hpp.
193template<class TypeTag>
194template <class Solver>
197step(const SimulatorTimer& simulator_timer,
198 Solver& solver,
199 const bool is_event,
200 const TuningUpdateCallback& tuning_updater)
201{
202 SubStepper<Solver> sub_stepper{
203 *this, simulator_timer, solver, is_event, tuning_updater,
204 };
205 return sub_stepper.run();
206}
207
208template<class TypeTag>
209template<class Serializer>
210void
212serializeOp(Serializer& serializer)
213{
214 serializer(this->time_step_control_type_);
215 switch (this->time_step_control_type_) {
217 allocAndSerialize<HardcodedTimeStepControl>(serializer);
218 break;
220 allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
221 break;
223 allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
224 break;
226 allocAndSerialize<PIDTimeStepControl>(serializer);
227 break;
229 allocAndSerialize<General3rdOrderController>(serializer);
230 break;
231 }
232 serializer(this->restart_factor_);
233 serializer(this->growth_factor_);
234 serializer(this->max_growth_);
235 serializer(this->max_time_step_);
236 serializer(this->min_time_step_);
237 serializer(this->ignore_convergence_failure_);
238 serializer(this->solver_restart_max_);
239 serializer(this->solver_verbose_);
240 serializer(this->timestep_verbose_);
241 serializer(this->suggested_next_timestep_);
242 serializer(this->full_timestep_initially_);
243 serializer(this->timestep_after_event_);
244 serializer(this->use_newton_iteration_);
245 serializer(this->min_time_step_before_shutting_problematic_wells_);
246}
247
248template<class TypeTag>
251report()
252{
253 return report_;
254}
255
256template<class TypeTag>
260{
261 return serializationTestObject_<HardcodedTimeStepControl>();
262}
263
264template<class TypeTag>
268{
269 return serializationTestObject_<PIDTimeStepControl>();
270}
271
272template<class TypeTag>
276{
277 return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
278}
279
280template<class TypeTag>
284{
285 return serializationTestObject_<SimpleIterationCountTimeStepControl>();
286}
287
288template<class TypeTag>
292{
293 return serializationTestObject_<General3rdOrderController>();
294}
295
296
297template<class TypeTag>
298void
300setSuggestedNextStep(const double x)
301{
302 this->suggested_next_timestep_ = x;
303}
304
305template<class TypeTag>
306double
308suggestedNextStep() const
309{
310 return this->suggested_next_timestep_;
311}
312
313template<class TypeTag>
316timeStepControl() const
317{
318 return *this->time_step_control_;
319}
320
321
322template<class TypeTag>
323void
325updateNEXTSTEP(double max_next_tstep)
326{
327 // \Note Only update next suggested step if TSINIT was explicitly
328 // set in TUNING or NEXTSTEP is active.
329 if (max_next_tstep > 0) {
330 this->suggested_next_timestep_ = max_next_tstep;
331 }
332}
333
334template<class TypeTag>
335void
337updateTUNING(double max_next_tstep, const Tuning& tuning)
338{
339 this->restart_factor_ = tuning.TSFCNV;
340 this->growth_factor_ = tuning.TFDIFF;
341 this->max_growth_ = tuning.TSFMAX;
342 this->max_time_step_ = tuning.TSMAXZ;
343 updateNEXTSTEP(max_next_tstep);
344 this->timestep_after_event_ = tuning.TMAXWC;
345}
346
347/*********************************************
348 * Private methods of AdaptiveTimeStepping
349 * ******************************************/
350
351template<class TypeTag>
352template<class T, class Serializer>
353void
355allocAndSerialize(Serializer& serializer)
356{
357 if (!serializer.isSerializing()) {
358 this->time_step_control_ = std::make_unique<T>();
359 }
360 serializer(*static_cast<T*>(this->time_step_control_.get()));
361}
362
363template<class TypeTag>
364template<class T>
365bool
366AdaptiveTimeStepping<TypeTag>::
367castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const
368{
369 const T* lhs = static_cast<const T*>(this->time_step_control_.get());
370 const T* rhs = static_cast<const T*>(Rhs.time_step_control_.get());
371 return *lhs == *rhs;
372}
373
374template<class TypeTag>
375void
376AdaptiveTimeStepping<TypeTag>::
377maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step,
378 bool is_event)
379{
380 // init last time step as a fraction of the given time step
381 if (this->suggested_next_timestep_ < 0) {
382 this->suggested_next_timestep_ = this->restart_factor_ * original_time_step;
383 }
384
385 if (this->full_timestep_initially_) {
386 this->suggested_next_timestep_ = original_time_step;
387 }
388
389 // use seperate time step after event
390 if (is_event && this->timestep_after_event_ > 0) {
391 this->suggested_next_timestep_ = this->timestep_after_event_;
392 }
393}
394
395template<class TypeTag>
396template<class Controller>
397AdaptiveTimeStepping<TypeTag>
398AdaptiveTimeStepping<TypeTag>::
399serializationTestObject_()
400{
401 AdaptiveTimeStepping<TypeTag> result;
402
403 result.restart_factor_ = 1.0;
404 result.growth_factor_ = 2.0;
405 result.max_growth_ = 3.0;
406 result.max_time_step_ = 4.0;
407 result.min_time_step_ = 5.0;
408 result.ignore_convergence_failure_ = true;
409 result.solver_restart_max_ = 6;
410 result.solver_verbose_ = true;
411 result.timestep_verbose_ = true;
412 result.suggested_next_timestep_ = 7.0;
413 result.full_timestep_initially_ = true;
414 result.use_newton_iteration_ = true;
415 result.min_time_step_before_shutting_problematic_wells_ = 9.0;
416 result.time_step_control_type_ = Controller::Type;
417 result.time_step_control_ =
418 std::make_unique<Controller>(Controller::serializationTestObject());
419
420 return result;
421}
422
423/*********************************************
424 * Protected methods of AdaptiveTimeStepping
425 * ******************************************/
426
427template<class TypeTag>
429init_(const UnitSystem& unitSystem)
430{
431 std::tie(time_step_control_type_,
432 time_step_control_,
433 use_newton_iteration_) = detail::createController(unitSystem);
434 // make sure growth factor is something reasonable
435 if (this->growth_factor_ < 1.0) {
436 OPM_THROW(std::runtime_error,
437 "Growth factor cannot be less than 1.");
438 }
439}
440
441
442
443/************************************************
444 * Private class SubStepper public methods
445 ************************************************/
446
447template<class TypeTag>
448template<class Solver>
450SubStepper(AdaptiveTimeStepping<TypeTag>& adaptive_time_stepping,
451 const SimulatorTimer& simulator_timer,
452 Solver& solver,
453 const bool is_event,
454 const TuningUpdateCallback& tuning_updater)
455 : adaptive_time_stepping_{adaptive_time_stepping}
456 , simulator_timer_{simulator_timer}
457 , solver_{solver}
458 , is_event_{is_event}
459 , tuning_updater_{tuning_updater}
460{
461}
462
463template<class TypeTag>
464template<class Solver>
465AdaptiveTimeStepping<TypeTag>&
466AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
467getAdaptiveTimerStepper()
468{
469 return adaptive_time_stepping_;
470}
471
472template<class TypeTag>
473template<class Solver>
474SimulatorReport
475AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
476run()
477{
478#ifdef RESERVOIR_COUPLING_ENABLED
479 if (isReservoirCouplingSlave_() && reservoirCouplingSlave_().activated()) {
480 return runStepReservoirCouplingSlave_();
481 }
482 else if (isReservoirCouplingMaster_() && reservoirCouplingMaster_().activated()) {
483 return runStepReservoirCouplingMaster_();
484 }
485 else {
486 return runStepOriginal_();
487 }
488#else
489 return runStepOriginal_();
490#endif
491}
492
493/************************************************
494 * Private class SubStepper private methods
495 ************************************************/
496
497#ifdef RESERVOIR_COUPLING_ENABLED
498template<class TypeTag>
499template<class Solver>
500bool
501AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
502isReservoirCouplingMaster_() const
503{
504 return this->solver_.model().simulator().reservoirCouplingMaster() != nullptr;
505}
506
507template<class TypeTag>
508template<class Solver>
509bool
510AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
511isReservoirCouplingSlave_() const
512{
513 return this->solver_.model().simulator().reservoirCouplingSlave() != nullptr;
514}
515#endif
516
517template<class TypeTag>
518template<class Solver>
519void
520AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
521maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step)
522{
523 this->adaptive_time_stepping_.maybeModifySuggestedTimeStepAtBeginningOfReportStep_(
524 original_time_step, this->is_event_
525 );
526
527 if (this->adaptive_time_stepping_.time_step_control_type_ != TimeStepControlType::HardCodedTimeStep) {
528 return;
529 }
530
531 struct ZeroRelativeChange final : RelativeChangeInterface {
532 double relativeChange() const override { return 0.0; }
533 } zero_relative_change;
534
535 const auto* hardcoded_control = static_cast<const HardcodedTimeStepControl*>(
536 this->adaptive_time_stepping_.time_step_control_.get());
537 AdaptiveSimulatorTimer report_step_timer{
538 this->simulator_timer_.startDateTime(),
539 original_time_step,
540 this->simulator_timer_.simulationTimeElapsed(),
541 original_time_step,
542 this->simulator_timer_.reportStepNum(),
543 maxTimeStep_()
544 };
545
546 const double hardcoded_initial_step = hardcoded_control->computeTimeStepSize(
547 original_time_step,
548 0,
549 zero_relative_change,
550 report_step_timer);
551 if (std::isfinite(hardcoded_initial_step) && hardcoded_initial_step > 0.0) {
552 this->adaptive_time_stepping_.setSuggestedNextStep(hardcoded_initial_step);
553 }
554}
555
556// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicit::runStep()
557// It has to be called for each substep since TUNING might have been changed for next sub step due
558// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
559template<class TypeTag>
560template<class Solver>
561bool
562AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
563maybeUpdateTuning_(double elapsed, double substep_length, int sub_step_number) const
564{
565 return this->tuning_updater_(elapsed, substep_length, sub_step_number);
566}
567
568template<class TypeTag>
569template<class Solver>
570double
571AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
572maxTimeStep_() const
573{
574 return this->adaptive_time_stepping_.max_time_step_;
575}
576
577template <class TypeTag>
578template <class Solver>
579SimulatorReport
580AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
581runStepOriginal_()
582{
583 const auto elapsed = this->simulator_timer_.simulationTimeElapsed();
584 const auto original_time_step = this->simulator_timer_.currentStepLength();
585 const auto report_step = this->simulator_timer_.reportStepNum();
586 maybeUpdateTuning_(elapsed, suggestedNextTimestep_(), /*substep=*/0);
587 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(original_time_step);
588
589 AdaptiveSimulatorTimer substep_timer{
590 this->simulator_timer_.startDateTime(),
591 original_time_step,
592 elapsed,
593 suggestedNextTimestep_(),
594 report_step,
595 maxTimeStep_()
596 };
597 SubStepIteration<Solver> substepIteration{*this, substep_timer, original_time_step, /*final_step=*/true};
598 return substepIteration.run();
599}
600
601template <class TypeTag>
602template <class Solver>
603double
604AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
605suggestedNextTimestep_() const
606{
607 return this->adaptive_time_stepping_.suggestedNextStep();
608}
609
610
611#ifdef RESERVOIR_COUPLING_ENABLED
612// Throw if the slave has already been terminated by the master. A terminated slave has
613// disconnected its intercommunicator, so running another coupled substep loop would issue
614// an MPI_Recv on a null communicator and abort the job. The run loop in
615// SimulatorFullyImplicit::runStep() stops before this can happen, so this guards
616// against future regressions of that logic.
617template <class TypeTag>
618template <class Solver>
619void
620AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
621checkIfSlaveIsTerminated_()
622{
623 if (reservoirCouplingSlave_().terminated()) {
624 OPM_THROW(ReservoirCouplingError,
625 "Internal error: attempt to run a coupled substep loop after the slave "
626 "has been terminated by the master process");
627 }
628}
629
630// Pick the master's sync-step length for the next outer-loop iteration of
631// `runStepReservoirCouplingMaster_()`. Includes the chop against slave-
632// report dates and emits the user-visible log line. See the block comment
633// above `runStepReservoirCouplingMaster_()` for TSYNC/RSYNC mode definitions.
634//
635// TSYNC (default):
636// Each iteration is one master time step. Start from the master's
637// adaptive suggestion, capped at the max time step and at the time
638// remaining until the report-step end, then chop.
639//
640// RSYNC:
641// The outer loop iterates per chunk between slave-report boundaries.
642// Re-use the previous iteration's chopped value as the candidate for
643// this iteration (the caller seeds the first iteration with
644// `original_time_step`), then chop.
645template <class TypeTag>
646template <class Solver>
647double
648AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
649getRcMasterSyncStepLength_(double prev_step,
650 double current_time,
651 double step_end_time)
652{
653 const bool sync_at_report_steps = reservoirCouplingMaster_().syncAtReportSteps();
654 double current_step_length;
655 if (sync_at_report_steps) {
656 current_step_length = prev_step;
657 } else {
658 const double remaining = step_end_time - current_time;
659 current_step_length = std::min({suggestedNextTimestep_(), maxTimeStep_(), remaining});
660 // NOTE: The substep timer is later constructed (in the caller) with
661 // current_step_length as its span, and after construction,
662 // substep_timer.currentStepLength() should return the same as
663 // current_step_length. The timer's constructor calls
664 // provideTimeStepEstimate() with suggestedNextTimestep_() as the dt_
665 // estimate and current_step_length as the span. Since we took
666 // current_step_length = min(suggestedNextTimestep_(), maxTimeStep_(), remaining)
667 // here (and only shrink it further via maybeChopSubStep below), the
668 // estimate is always >= the span, so the timer's snap branch (see
669 // AdaptiveSimulatorTimer::provideTimeStepEstimate) fires and clamps
670 // dt_ to current_step_length.
671 }
672 current_step_length = reservoirCouplingMaster_().maybeChopSubStep(current_step_length, current_time);
673 auto num_active = reservoirCouplingMaster_().numActivatedSlaves();
674 OpmLog::info(fmt::format(
675 "\nChoosing next sync time{} between master and {} active slave {}: {:.2f} days",
676 sync_at_report_steps ? " (RSYNC)" : "",
677 num_active, (num_active == 1 ? "process" : "processes"),
678 current_step_length / unit::day
679 ));
680 return current_step_length;
681}
682
683template <class TypeTag>
684template <class Solver>
685ReservoirCouplingMaster<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
686AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
687reservoirCouplingMaster_()
688{
689 return *(this->solver_.model().simulator().reservoirCouplingMaster());
690}
691
692template <class TypeTag>
693template <class Solver>
694ReservoirCouplingSlave<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
695AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
696reservoirCouplingSlave_()
697{
698 return *(this->solver_.model().simulator().reservoirCouplingSlave());
699}
700
701// Description of the reservoir coupling master and slave substep loop
702// -------------------------------------------------------------------
703// The master and slave processes attempt to reach the end of the report step using a series of substeps
704// (also called timesteps). Each substep has an upper limit roughly determined by a combination of the
705// keywords TUNING (through TSINIT, TSMAXZ), NEXSTEP, WCYCLE, and the start of the next report step
706// (the last substep has to coincide with this time). Note that NEXTSTEP can be updated from an
707// ACTIONX keyword. Although this comment focuses on the maximum substep limit, there is also a lower
708// limit on the substep length, and the substep sizes are adjusted automatically (or retried) based on
709// the convergence behavior of the solver and other criteria.
710//
711// The master chooses a "sync step" — the span between successive master<->slave exchanges — via one
712// of two modes, selected by the `--rescoup-sync-at-report-steps` CLI flag:
713//
714// TSYNC (default, flag=false)
715// Each master time step is a sync step. Every outer iteration of the master substep
716// loop redistributes group-rate targets to the slaves. This gives fine sync granularity and
717// matches the reference simulator's coupling semantics.
718//
719// RSYNC (flag=true)
720// Each chunk between slave-report-step boundaries is one sync step. The master's own adaptive
721// timestepping runs inside the chunk, but slaves only hear back from the master once per chunk.
722// Provided as a developer escape hatch.
723//
724// In both modes the sync step is also limited so as not to overshoot the next slave report date
725// (`maybeChopSubStep`). A second, drift-based limit (GRUPMAST item 4 / RCMASTS) is parsed but not
726// yet honored — tracked as a follow-up PR.
727//
728// Per sync step:
729// - Master receives each activated slave's next report date (`receiveNextReportDateFromSlaves`).
730// - Master picks the sync-step length via `getRcMasterSyncStepLength_()` (mode-dependent) and
731// chops it so as not to overshoot any slave report date.
732// - Master sends the chosen sync-step length to each slave (`sendNextTimeStepToSlaves`).
733//
734// The slaves use the sync-step end as a fixed point — a "mini report step" — that they reach via one
735// or more of their own substeps. If the slave's current report step extends beyond the received
736// sync step, the slave loops waiting for the master to send further sync steps; the loop ends when
737// the received sync step coincides with the end of the slave's current report step.
738
739// Master-side rescoup outer loop. See the block comment above for TSYNC
740// and RSYNC mode definitions.
741//
742// TSYNC (default):
743// Each outer iteration is one master time step. The sync-step
744// length is chosen fresh each iteration by `getRcMasterSyncStepLength_()`
745// from the master's adaptive suggestion (capped at the max time step and
746// at the remaining report-step time), then chopped against slave-report
747// dates. `SubStepIteration::run()` is still used but normally handles
748// only solver-driven chops within a single master time step.
749//
750// RSYNC:
751// The outer loop iterates once per chunk between slave-report boundaries.
752// The sync-step length starts at the full report step and only shrinks
753// via `maybeChopSubStep`; the master's adaptive timestepping runs inside
754// `SubStepIteration` as solver-driven substeps within each chunk.
755template <class TypeTag>
756template <class Solver>
757SimulatorReport
758AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
759runStepReservoirCouplingMaster_()
760{
761 int iteration = 0;
762 const double original_time_step = this->simulator_timer_.currentStepLength();
763 double current_time{this->simulator_timer_.simulationTimeElapsed()};
764 double step_end_time = current_time + original_time_step;
765 const double report_step_start_time = current_time;
766 int report_step_substep_offset = 0;
767 // In RSYNC mode this variable persists across outer iterations and
768 // carries the previously-chopped sync span into the next iteration. In
769 // TSYNC mode it is overwritten each iteration by
770 // `getRcMasterSyncStepLength_()`; the initial value is only consumed by
771 // the first helper call in RSYNC mode (as the candidate span for that
772 // iteration).
773 auto current_step_length = original_time_step;
774 auto report_step_idx = this->simulator_timer_.currentStepNum();
775 if (report_step_idx == 0 && iteration == 0) {
776 reservoirCouplingMaster_().initTimeStepping();
777 }
778 SimulatorReport report;
779 // The master needs to know which slaves have activated before it can start the substep loop
780 reservoirCouplingMaster_().maybeReceiveActivationHandshakeFromSlaves(current_time);
781 while (true) {
782 reservoirCouplingMaster_().sendDontTerminateSignalToSlaves(); // Tell the slaves to keep running.
783 reservoirCouplingMaster_().receiveNextReportDateFromSlaves();
784 const bool start_of_report_step = (iteration == 0);
785 if (start_of_report_step) {
786 reservoirCouplingMaster_().initStartOfReportStep(report_step_idx);
787 maybeUpdateTuning_(current_time, suggestedNextTimestep_(), /*substep=*/0);
788 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(original_time_step);
789 }
790 current_step_length = getRcMasterSyncStepLength_(
791 current_step_length, current_time, step_end_time);
792 reservoirCouplingMaster_().sendNextTimeStepToSlaves(current_step_length);
793 AdaptiveSimulatorTimer substep_timer{
794 this->simulator_timer_.startDateTime(),
795 /*stepLength=*/current_step_length,
796 /*elapsedTime=*/current_time,
797 /*timeStepEstimate=*/suggestedNextTimestep_(),
798 /*reportStep=*/this->simulator_timer_.reportStepNum(),
799 maxTimeStep_()
800 };
801 // Make the per-chunk timer log the enclosing report step's span and a
802 // cumulative substep counter (see AdaptiveSimulatorTimer "report step
803 // view"). Without this, log lines would show one sync chunk only.
804 substep_timer.setReportStepStartTime(report_step_start_time);
805 substep_timer.setReportStepTotalTime(step_end_time);
806 substep_timer.setReportStepSubstepOffset(report_step_substep_offset);
807 const bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq(
808 current_time + current_step_length, step_end_time
809 );
810 // Mark this as the first substep of the "sync" timestep. This flag controls
811 // whether master-slave data exchange should occur in beginTimeStep() in the well model.
812 // It will be cleared after the first runSubStep_() call.
813 reservoirCouplingMaster_().setFirstSubstepOfSyncTimestep(true);
814 // After the first master substep completes, timeStepSucceeded() will
815 // block until slaves finish the sync step and send production data.
816 // This ensures correct summary output for all subsequent substeps.
817 reservoirCouplingMaster_().setNeedsSlaveDataReceive(true);
818 SubStepIteration<Solver> substepIteration{*this, substep_timer, current_step_length, final_step};
819 const auto sub_steps_report = substepIteration.run();
820 report += sub_steps_report;
821 report_step_substep_offset += substep_timer.currentStepNum();
822 current_time += current_step_length;
823 if (final_step) {
824 break;
825 }
826 iteration++;
827 }
828 return report;
829}
830
831template <class TypeTag>
832template <class Solver>
833SimulatorReport
834AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
835runStepReservoirCouplingSlave_()
836{
837 checkIfSlaveIsTerminated_();
838 int iteration = 0;
839 const double original_time_step = this->simulator_timer_.currentStepLength();
840 double current_time{this->simulator_timer_.simulationTimeElapsed()};
841 double step_end_time = current_time + original_time_step;
842 const double report_step_start_time = current_time;
843 int report_step_substep_offset = 0;
844 SimulatorReport report;
845 auto report_step_idx = this->simulator_timer_.currentStepNum();
846 if (report_step_idx == 0 && iteration == 0) {
847 reservoirCouplingSlave_().initTimeStepping();
848 }
849 while (true) {
850 bool start_of_report_step = (iteration == 0);
851 if (reservoirCouplingSlave_().maybeReceiveTerminateSignalFromMaster()) {
852 // Call MPI_Comm_disconnect() to terminate the MPI communicator, etc..
853 break;
854 }
855 reservoirCouplingSlave_().sendNextReportDateToMasterProcess();
856 const auto timestep = reservoirCouplingSlave_().receiveNextTimeStepFromMaster();
857 if (start_of_report_step) {
858 maybeUpdateTuning_(current_time, suggestedNextTimestep_(), /*substep=*/0);
859 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(timestep);
860 }
861 AdaptiveSimulatorTimer substep_timer{
862 this->simulator_timer_.startDateTime(),
863 /*step_length=*/timestep,
864 /*elapsed_time=*/current_time,
865 /*time_step_estimate=*/suggestedNextTimestep_(),
866 this->simulator_timer_.reportStepNum(),
867 maxTimeStep_()
868 };
869 // Make the per-chunk timer log the enclosing report step's span and a
870 // cumulative substep counter (see AdaptiveSimulatorTimer "report step
871 // view"). Without this, log lines would show one sync chunk only.
872 substep_timer.setReportStepStartTime(report_step_start_time);
873 substep_timer.setReportStepTotalTime(step_end_time);
874 substep_timer.setReportStepSubstepOffset(report_step_substep_offset);
875 const bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq(
876 current_time + timestep, step_end_time
877 );
878 // Mark this as the first substep of the "sync" timestep. This flag controls
879 // whether master-slave data exchange should occur in beginTimeStep() in the well model.
880 // It will be cleared after the first runSubStep_() call.
881 reservoirCouplingSlave_().setFirstSubstepOfSyncTimestep(true);
882 SubStepIteration<Solver> substepIteration{*this, substep_timer, timestep, final_step};
883 const auto sub_steps_report = substepIteration.run();
884 report += sub_steps_report;
885 report_step_substep_offset += substep_timer.currentStepNum();
886 current_time += timestep;
887 if (final_step) {
888 break;
889 }
890 iteration++;
891 }
892 return report;
893}
894#endif // RESERVOIR_COUPLING_ENABLED
895
896
897/************************************************
898 * Private class SubStepIteration public methods
899 ************************************************/
900
901template<class TypeTag>
902template<class Solver>
903AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
904SubStepIteration(SubStepper<Solver>& substepper,
905 AdaptiveSimulatorTimer& substep_timer,
906 const double original_time_step,
907 bool final_step)
908 : substepper_{substepper}
909 , substep_timer_{substep_timer}
910 , original_time_step_{original_time_step}
911 , final_step_{final_step}
912 , adaptive_time_stepping_{substepper.getAdaptiveTimerStepper()}
913{
914}
915
916template <class TypeTag>
917template <class Solver>
918SimulatorReport
919AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
920run()
921{
922 auto& simulator = solver_().model().simulator();
923 auto& problem = simulator.problem();
924 // counter for solver restarts
925 int restarts = 0;
926 SimulatorReport report;
927
928 // sub step time loop
929 while (!this->substep_timer_.done()) {
930 // if we just chopped the timestep due to convergence i.e. restarts>0
931 // we dont want to update the next timestep based on Tuning
932 if (restarts == 0) {
933 maybeUpdateTuningAndTimeStep_();
934 }
935 const double dt = this->substep_timer_.currentStepLength();
936 if (timeStepVerbose_()) {
937 detail::logTimer(this->substep_timer_);
938 }
939
940 maybeUpdateLastSubstepOfSyncTimestep_(dt); // Needed for reservoir coupling
941 auto substep_report = runSubStep_();
942 markFirstSubStepAsFinished_(); // Needed for reservoir coupling
943
944 if (substep_report.converged || checkContinueOnUnconvergedSolution_(dt)) {
945 Dune::Timer perfTimer;
946 perfTimer.start();
947 // Pass substep to eclwriter for summary output
948 problem.setSubStepReport(substep_report);
949 auto& full_report = adaptive_time_stepping_.report();
950 full_report += substep_report;
951 problem.setSimulationReport(full_report);
952 problem.endTimeStep();
953 substep_report.pre_post_time += perfTimer.stop();
954
955 report += substep_report;
956
957 OPM_TIMEBLOCK(convergenceSucceeded);
958 ++this->substep_timer_; // advance by current dt
959
960 const int iterations = getNumIterations_(substep_report);
961 auto dt_estimate = timeStepControlComputeEstimate_(
962 dt, iterations, this->substep_timer_);
963
964 assert(dt_estimate > 0);
965 dt_estimate = maybeRestrictTimeStepGrowth_(dt, dt_estimate, restarts);
966 restarts = 0; // solver converged, reset restarts counter
967
968 maybeReportSubStep_(substep_report);
969 if (this->final_step_ && this->substep_timer_.done()) {
970 // if the time step is done we do not need to write it as this will be done
971 // by the simulator anyway.
972 }
973 else {
974 report.success.output_write_time += writeOutput_();
975 }
976
977 // set new time step length
978 checkTimeStepCanAdvance_(dt, dt_estimate);
979 setTimeStep_(dt_estimate);
980
981 report.success.converged = this->substep_timer_.done();
982 this->substep_timer_.setLastStepFailed(false);
983 }
984 else { // in case of no convergence or time step tolerance test failure
985 OPM_TIMEBLOCK(convergenceFailed);
986 report += substep_report;
987 this->substep_timer_.setLastStepFailed(true);
988 checkTimeStepMaxRestartLimit_(restarts);
989
990 double new_time_step = restartFactor_() * dt;
991 if (substep_report.time_step_rejected) {
992 const double tol = Parameters::Get<Parameters::TimeStepControlTolerance>();
993 const double safetyFactor = Parameters::Get<Parameters::TimeStepControlSafetyFactor>();
994 const double temp_time_step = std::sqrt(safetyFactor * tol / solver_().model().relativeChange()) * dt;
995 if (temp_time_step < dt) { // added in case suggested time step is not a reduction
996 new_time_step = temp_time_step;
997 }
998 }
999 checkTimeStepMinLimit_(new_time_step);
1000 bool wells_shut = false;
1001 if (new_time_step > minTimeStepBeforeClosingWells_()) {
1002 chopTimeStep_(new_time_step);
1003 } else {
1004 wells_shut = chopTimeStepOrCloseFailingWells_(new_time_step);
1005 }
1006 if (wells_shut) {
1007 setTimeStep_(dt); // retry the old timestep
1008 }
1009 else {
1010 restarts++; // only increase if no wells were shut
1011 }
1012 }
1013 problem.setNextTimeStepSize(this->substep_timer_.currentStepLength());
1014 }
1015 updateSuggestedNextStep_();
1016 return report;
1017}
1018
1019
1020/************************************************
1021 * Private class SubStepIteration private methods
1022 ************************************************/
1023
1024
1025template<class TypeTag>
1026template<class Solver>
1027bool
1028AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1029checkContinueOnUnconvergedSolution_(double dt) const
1030{
1031 const bool continue_on_uncoverged_solution = ignoreConvergenceFailure_() && dt <= minTimeStep_();
1032 if (continue_on_uncoverged_solution && solverVerbose_()) {
1033 // NOTE: This method is only called if the solver failed to converge.
1034 const auto msg = fmt::format(
1035 "Solver failed to converge but timestep {} is smaller or equal to {}\n"
1036 "which is the minimum threshold given by option --solver-min-time-step\n",
1037 dt, minTimeStep_()
1038 );
1039 OpmLog::problem(msg);
1040 }
1041 return continue_on_uncoverged_solution;
1042}
1043
1044template<class TypeTag>
1045template<class Solver>
1046void
1047AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1048checkTimeStepMaxRestartLimit_(const int restarts) const
1049{
1050 // If we have restarted (i.e. cut the timestep) too
1051 // many times, we have failed and throw an exception.
1052 if (restarts >= solverRestartMax_()) {
1053 const auto msg = fmt::format(
1054 fmt::runtime("Solver failed to converge after cutting timestep {} times."), restarts
1055 );
1056 if (solverVerbose_()) {
1057 OpmLog::error(msg);
1058 }
1059 // Use throw directly to prevent file and line
1060 throw TimeSteppingBreakdown{msg};
1061 }
1062}
1063
1064template<class TypeTag>
1065template<class Solver>
1066void
1067AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1068checkTimeStepCanAdvance_(const double current_time_step,
1069 const double new_time_step) const
1070{
1071 // Both dt floors - the minimum step size and the restart count - sit on
1072 // the failure path. A run whose substeps all converge, but for which the
1073 // controller keeps proposing ever smaller steps, consults neither: the
1074 // step size shrinks without bound while the elapsed time stands still,
1075 // and the run never ends. Hold the accepted path to the same minimum.
1076 if (new_time_step >= minTimeStep_()) {
1077 return;
1078 }
1079
1080 // A sub-minimum step is valid when the report step is already complete, or
1081 // when timestep control is increasing it. Reject only a shrinking
1082 // sub-minimum step, which is the one that stalls short of the report-step
1083 // boundary.
1084 if (this->substep_timer_.done() || new_time_step >= current_time_step) {
1085 return;
1086 }
1087
1088 const auto msg =
1089 fmt::format("Time step control proposed a step of {:.3E} DAYS, below the "
1090 "minimum of {:.3E} DAYS, while every substep converges. The "
1091 "run is not advancing past {:.6E} DAYS.",
1092 new_time_step / 86400.0,
1093 minTimeStep_() / 86400.0,
1094 this->substep_timer_.simulationTimeElapsed() / 86400.0);
1095
1096 if (solverVerbose_()) {
1097 OpmLog::error(msg);
1098 }
1099
1100 // Use throw directly to prevent file and line
1101 throw TimeSteppingBreakdown{msg};
1102}
1103
1104template<class TypeTag>
1105template<class Solver>
1106void
1107AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1108checkTimeStepMinLimit_(const double new_time_step) const
1109{
1110 using Meas = UnitSystem::measure;
1111 // If we have restarted (i.e. cut the timestep) too
1112 // much, we have failed and throw an exception.
1113 if (new_time_step < minTimeStep_()) {
1114 std::string msg = "Solver failed to converge after cutting timestep to ";
1115 if (Parameters::Get<Parameters::EnableTuning>()) {
1116 const UnitSystem& unit_system = solver_().model().simulator().vanguard().eclState().getDeckUnitSystem();
1117 msg += fmt::format(
1118 "{:.3E} {}\nwhich is the minimum threshold given by the TUNING keyword\n",
1119 unit_system.from_si(Meas::time, minTimeStep_()),
1120 unit_system.name(Meas::time)
1121 );
1122 }
1123 else {
1124 msg += fmt::format(
1125 "{:.3E} DAYS\nwhich is the minimum threshold given by option --solver-min-time-step\n",
1126 minTimeStep_() / 86400.0
1127 );
1128 }
1129 if (solverVerbose_()) {
1130 OpmLog::error(msg);
1131 }
1132 // Use throw directly to prevent file and line
1133 throw TimeSteppingBreakdown{msg};
1134 }
1135}
1136
1137template<class TypeTag>
1138template<class Solver>
1139void
1140AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1141chopTimeStep_(const double new_time_step)
1142{
1143 setTimeStep_(new_time_step);
1144 if (solverVerbose_()) {
1145 const auto msg = fmt::format(fmt::runtime("{}\nTimestep chopped to {} days\n"),
1146 this->cause_of_failure_,
1147 unit::convert::to(this->substep_timer_.currentStepLength(), unit::day));
1148 OpmLog::problem(msg);
1149 }
1150}
1151
1152template<class TypeTag>
1153template<class Solver>
1154bool
1155AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1156chopTimeStepOrCloseFailingWells_(const double new_time_step)
1157{
1158 bool wells_shut = false;
1159 // We are below the threshold, and will check if there are any
1160 // wells that fails repeatedly (that means that it fails in the last three steps)
1161 // we should close rather than chopping again.
1162 // If we already have chopped the timestep two times that is
1163 // new_time_step < minTimeStepBeforeClosingWells_()*restartFactor_()*restartFactor_()
1164 // We also shut wells that fails only on this step.
1165 const bool requireRepeatedFailures =
1166 new_time_step > (minTimeStepBeforeClosingWells_() * restartFactor_() * restartFactor_());
1167 const std::set<std::string> failing_wells =
1168 detail::consistentlyFailingWells(solver_().model().stepReports(), requireRepeatedFailures);
1169
1170 if (failing_wells.empty()) {
1171 // Found no wells to close, chop the timestep
1172 chopTimeStep_(new_time_step);
1173 } else {
1174 // Close all consistently failing wells that are not under group control
1175 std::vector<std::string> shut_wells;
1176 for (const auto& well : failing_wells) {
1177 const bool was_shut =
1178 solver_().model().wellModel().forceShutWellByName(well,
1179 this->substep_timer_.simulationTimeElapsed(),
1180 /*dont_shut_grup_wells =*/ true);
1181 if (was_shut) {
1182 shut_wells.push_back(well);
1183 }
1184 }
1185 // If no wells are closed we also try to shut wells under group control
1186 if (shut_wells.empty()) {
1187 for (const auto& well : failing_wells) {
1188 const bool was_shut =
1189 solver_().model().wellModel().forceShutWellByName(well,
1190 this->substep_timer_.simulationTimeElapsed(),
1191 /*dont_shut_grup_wells =*/ false);
1192 if (was_shut) {
1193 shut_wells.push_back(well);
1194 }
1195 }
1196 }
1197 // If still no wells are closed we must fall back to chopping again
1198 if (shut_wells.empty()) {
1199 chopTimeStep_(new_time_step);
1200 } else {
1201 wells_shut = true;
1202 if (solverVerbose_()) {
1203 const std::string msg =
1204 fmt::format(fmt::runtime("\nProblematic well(s) were shut: {}"
1205 "(retrying timestep)\n"),
1206 fmt::join(shut_wells, " "));
1207 OpmLog::problem(msg);
1208 }
1209 }
1210 }
1211 return wells_shut;
1212}
1213
1214template<class TypeTag>
1215template<class Solver>
1216boost::posix_time::ptime
1217AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1218currentDateTime_() const
1219{
1220 return simulatorTimer_().currentDateTime();
1221}
1222
1223template<class TypeTag>
1224template<class Solver>
1225int
1226AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1227getNumIterations_(const SimulatorReportSingle &substep_report) const
1228{
1229 if (useNewtonIteration_()) {
1230 return substep_report.total_newton_iterations;
1231 }
1232 else {
1233 return substep_report.total_linear_iterations;
1234 }
1235}
1236
1237template<class TypeTag>
1238template<class Solver>
1239double
1240AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1241growthFactor_() const
1242{
1243 return this->adaptive_time_stepping_.growth_factor_;
1244}
1245
1246template<class TypeTag>
1247template<class Solver>
1248bool
1249AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1250ignoreConvergenceFailure_() const
1251{
1252 return adaptive_time_stepping_.ignore_convergence_failure_;
1253}
1254
1255template<class TypeTag>
1256template<class Solver>
1257bool
1258AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1259isReservoirCouplingMaster_() const
1260{
1261 return this->substepper_.isReservoirCouplingMaster_();
1262}
1263
1264template<class TypeTag>
1265template<class Solver>
1266bool
1267AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1268isReservoirCouplingSlave_() const
1269{
1270 return this->substepper_.isReservoirCouplingSlave_();
1271}
1272
1273template<class TypeTag>
1274template<class Solver>
1275void
1276AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1277markFirstSubStepAsFinished_() const
1278{
1279#ifdef RESERVOIR_COUPLING_ENABLED
1280 // Clear the first-substep flag after the first runSubStep_() call.
1281 // This ensures that master-slave synchronization only happens once per sync timestep,
1282 // not on retry attempts after convergence-driven timestep chops.
1283 if (isReservoirCouplingMaster_()) {
1284 reservoirCouplingMaster_().setFirstSubstepOfSyncTimestep(false);
1285 }
1286 else if (isReservoirCouplingSlave_()) {
1287 reservoirCouplingSlave_().setFirstSubstepOfSyncTimestep(false);
1288 }
1289#endif
1290 return;
1291}
1292
1293template<class TypeTag>
1294template<class Solver>
1295double
1296AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1297maxGrowth_() const
1298{
1299 return this->adaptive_time_stepping_.max_growth_;
1300}
1301
1302template<class TypeTag>
1303template<class Solver>
1304void
1305AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1306maybeReportSubStep_(SimulatorReportSingle substep_report) const
1307{
1308 if (timeStepVerbose_()) {
1309 std::ostringstream ss;
1310 substep_report.reportStep(ss);
1311 OpmLog::info(ss.str());
1312 }
1313}
1314
1315template<class TypeTag>
1316template<class Solver>
1317double
1318AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1319maybeRestrictTimeStepGrowth_(const double dt, double dt_estimate, const int restarts) const
1320{
1321 // limit the growth of the timestep size by the growth factor
1322 dt_estimate = std::min(dt_estimate, double(maxGrowth_() * dt));
1323 assert(dt_estimate > 0);
1324 // further restrict time step size growth after convergence problems
1325 if (restarts > 0) {
1326 dt_estimate = std::min(growthFactor_() * dt, dt_estimate);
1327 }
1328
1329 return dt_estimate;
1330}
1331
1332
1333template<class TypeTag>
1334template<class Solver>
1335void
1336AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1337maybeUpdateLastSubstepOfSyncTimestep_([[maybe_unused]] const double dt)
1338{
1339#ifdef RESERVOIR_COUPLING_ENABLED
1340 // For reservoir coupling slaves: predict if this substep will complete
1341 // the sync timestep. If so, timeStepSucceeded() will send production
1342 // data to the master (which is blocking on receive after its first substep).
1343 // This is used for summary data synchronization between slaves and master.
1344 if (isReservoirCouplingSlave_()) {
1346 this->substep_timer_.simulationTimeElapsed() + dt,
1347 this->substep_timer_.totalTime()
1348 );
1349 reservoirCouplingSlave_().setLastSubstepOfSyncTimestep(is_last);
1350 }
1351#endif
1352}
1353
1354// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicit::runStep()
1355// It has to be called for each substep since TUNING might have been changed for next sub step due
1356// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
1357template<class TypeTag>
1358template<class Solver>
1359void
1360AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1361maybeUpdateTuningAndTimeStep_()
1362{
1363 // TODO: This function is currently only called if NEXTSTEP is activated from ACTIONX or
1364 // if the WCYCLE keyword needs to modify the current timestep. So this method should rather
1365 // be named maybeUpdateTimeStep_() or similar, since it should not update the tuning. However,
1366 // the current definition of the maybeUpdateTuning_() callback is actually calling
1367 // adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning) which is updating the tuning
1368 // see SimulatorFullyImplicit::runStep() for more details.
1369 const auto old_value = suggestedNextTimestep_();
1370 if (this->substepper_.maybeUpdateTuning_(this->substep_timer_.simulationTimeElapsed(),
1371 this->substep_timer_.currentStepLength(),
1372 this->substep_timer_.currentStepNum()))
1373 {
1374 // Either NEXTSTEP and WCYCLE wants to change the current time step, but they cannot
1375 // change the current time step directly. Instead, they change the suggested next time step
1376 // by calling updateNEXTSTEP() via the maybeUpdateTuning() callback. We now need to update
1377 // the current time step to the new suggested time step and reset the suggested time step
1378 // to the old value.
1379 setTimeStep_(suggestedNextTimestep_());
1380 setSuggestedNextStep_(old_value);
1381 }
1382}
1383
1384template<class TypeTag>
1385template<class Solver>
1386double
1387AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1388minTimeStepBeforeClosingWells_() const
1389{
1390 return this->adaptive_time_stepping_.min_time_step_before_shutting_problematic_wells_;
1391}
1392
1393template<class TypeTag>
1394template<class Solver>
1395double
1396AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1397minTimeStep_() const
1398{
1399 return this->adaptive_time_stepping_.min_time_step_;
1400}
1401
1402#ifdef RESERVOIR_COUPLING_ENABLED
1403template<class TypeTag>
1404template<class Solver>
1405ReservoirCouplingMaster<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
1406AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1407reservoirCouplingMaster_() const
1408{
1409 return this->substepper_.reservoirCouplingMaster_();
1410}
1411
1412template<class TypeTag>
1413template<class Solver>
1414ReservoirCouplingSlave<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
1415AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1416reservoirCouplingSlave_() const
1417{
1418 return this->substepper_.reservoirCouplingSlave_();
1419}
1420#endif
1421
1422template<class TypeTag>
1423template<class Solver>
1424double
1425AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1426restartFactor_() const
1427{
1428 return this->adaptive_time_stepping_.restart_factor_;
1429}
1430
1431template<class TypeTag>
1432template<class Solver>
1433SimulatorReportSingle
1434AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1435runSubStep_()
1436{
1437 OPM_TIMEFUNCTION();
1438 SimulatorReportSingle substep_report;
1439
1440 auto handleFailure = [this, &substep_report]
1441 (const std::string& failure_reason, const std::exception& e, bool log_exception = true)
1442 {
1443 substep_report = solver_().failureReport();
1444 // Keep the exception text: it carries the actual reason (e.g. which
1445 // equation had a too-large residual), which otherwise only surfaces
1446 // at debug verbosity while the PRT shows just the generic category.
1447 this->cause_of_failure_ = failure_reason;
1448 if (const std::string what = e.what();
1449 !what.empty() && what != failure_reason)
1450 {
1451 this->cause_of_failure_ += " (" + what + ")";
1452 }
1453 if (log_exception && solverVerbose_()) {
1454 OpmLog::debug(std::string("Caught Exception: ") + e.what());
1455 }
1456 };
1457
1458 try {
1459 substep_report = solver_().step(this->substep_timer_, &this->adaptive_time_stepping_.timeStepControl());
1460 if (solverVerbose_()) {
1461 // report number of linear iterations
1462 OpmLog::debug("Overall linear iterations used: "
1463 + std::to_string(substep_report.total_linear_iterations));
1464 }
1465 }
1466 catch (const TooManyIterations& e) {
1467 handleFailure("Solver convergence failure - Iteration limit reached", e);
1468 }
1469 catch (const TimeSteppingBreakdown& e) {
1470 handleFailure(e.what(), e);
1471 }
1472 catch (const ConvergenceMonitorFailure& e) {
1473 handleFailure("Convergence monitor failure", e, /*log_exception=*/false);
1474 }
1475 catch (const LinearSolverProblem& e) {
1476 handleFailure("Linear solver convergence failure", e);
1477 }
1478 catch (const NumericalProblem& e) {
1479 handleFailure("Solver convergence failure - Numerical problem encountered", e);
1480 }
1481 catch (const std::runtime_error& e) {
1482 handleFailure("Runtime error encountered", e);
1483 }
1484 catch (const Dune::ISTLError& e) {
1485 handleFailure("ISTL error - Time step too large", e);
1486 }
1487 catch (const Dune::MatrixBlockError& e) {
1488 handleFailure("Matrix block error", e);
1489 }
1490
1491 return substep_report;
1492}
1493
1494template<class TypeTag>
1495template<class Solver>
1496void
1497AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1498setTimeStep_(double dt_estimate)
1499{
1500 this->substep_timer_.provideTimeStepEstimate(dt_estimate);
1501}
1502
1503template<class TypeTag>
1504template<class Solver>
1505Solver&
1506AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1507solver_() const
1508{
1509 return this->substepper_.solver_;
1510}
1511
1512
1513template<class TypeTag>
1514template<class Solver>
1515int
1516AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1517solverRestartMax_() const
1518{
1519 return this->adaptive_time_stepping_.solver_restart_max_;
1520}
1521
1522template<class TypeTag>
1523template<class Solver>
1524void
1525AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1526setSuggestedNextStep_(double step)
1527{
1528 this->adaptive_time_stepping_.setSuggestedNextStep(step);
1529}
1530
1531template <class TypeTag>
1532template <class Solver>
1533const SimulatorTimer&
1534AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1535simulatorTimer_() const
1536{
1537 return this->substepper_.simulator_timer_;
1538}
1539
1540template <class TypeTag>
1541template <class Solver>
1542bool
1543AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1544solverVerbose_() const
1545{
1546 return this->adaptive_time_stepping_.solver_verbose_;
1547}
1548
1549template<class TypeTag>
1550template<class Solver>
1551boost::posix_time::ptime
1552AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1553startDateTime_() const
1554{
1555 return simulatorTimer_().startDateTime();
1556}
1557
1558template <class TypeTag>
1559template <class Solver>
1560double
1561AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1562suggestedNextTimestep_() const
1563{
1564 return this->adaptive_time_stepping_.suggestedNextStep();
1565}
1566
1567template <class TypeTag>
1568template <class Solver>
1569double
1570AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1571timeStepControlComputeEstimate_(const double dt, const int iterations,
1572 const AdaptiveSimulatorTimer& substepTimer) const
1573{
1574 // create object to compute the time error, simply forwards the call to the model
1575 const SolutionTimeErrorSolverWrapper<Solver> relative_change{solver_()};
1576 return this->adaptive_time_stepping_.time_step_control_->computeTimeStepSize(
1577 dt, iterations, relative_change, substepTimer);
1578}
1579
1580template <class TypeTag>
1581template <class Solver>
1582bool
1583AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1584timeStepVerbose_() const
1585{
1586 return this->adaptive_time_stepping_.timestep_verbose_;
1587}
1588
1589// The suggested time step is the stepsize that will be used as a first try for
1590// the next sub step. It is updated at the end of each substep. It can also be
1591// updated by the TUNING or NEXTSTEP keywords at the beginning of each report step or
1592// at the beginning of each substep by the ACTIONX keyword (via NEXTSTEP), this is
1593// done by the maybeUpdateTuning_() method which is called at the beginning of each substep
1594// (and the begginning of each report step). Note that the WCYCLE keyword can also update the
1595// suggested time step via the maybeUpdateTuning_() method.
1596template <class TypeTag>
1597template <class Solver>
1598void
1599AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1600updateSuggestedNextStep_()
1601{
1602 auto suggested_next_step = this->substep_timer_.currentStepLength();
1603 if (! std::isfinite(suggested_next_step)) { // check for NaN
1604 suggested_next_step = this->original_time_step_;
1605 }
1606 if (timeStepVerbose_()) {
1607 std::ostringstream ss;
1608 this->substep_timer_.report(ss);
1609 ss << "Suggested next step size = "
1610 << unit::convert::to(suggested_next_step, unit::day) << " (days)" << std::endl;
1611 OpmLog::debug(ss.str());
1612 }
1613 setSuggestedNextStep_(suggested_next_step);
1614}
1615
1616template <class TypeTag>
1617template <class Solver>
1618bool
1619AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1620useNewtonIteration_() const
1621{
1622 return this->adaptive_time_stepping_.use_newton_iteration_;
1623}
1624
1625template <class TypeTag>
1626template <class Solver>
1627double
1628AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1629writeOutput_() const
1630{
1631 time::StopWatch perf_timer;
1632 perf_timer.start();
1633 auto& problem = solver_().model().simulator().problem();
1634 problem.writeOutput(true);
1635 return perf_timer.secsSinceStart();
1636}
1637
1638/************************************************
1639 * Private class SolutionTimeErrorSolverWrapper
1640 * **********************************************/
1641
1642template<class TypeTag>
1643template<class Solver>
1644AdaptiveTimeStepping<TypeTag>::
1645SolutionTimeErrorSolverWrapper<Solver>::
1646SolutionTimeErrorSolverWrapper(const Solver& solver)
1647 : solver_{solver}
1648{}
1649
1650template<class TypeTag>
1651template<class Solver>
1652double AdaptiveTimeStepping<TypeTag>::SolutionTimeErrorSolverWrapper<Solver>::relativeChange() const
1653{
1654 // returns: || u^n+1 - u^n || / || u^n+1 ||
1655 return solver_.model().relativeChange();
1656}
1657
1658} // namespace Opm
1659
1660#endif // OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
Adaptive time-stepping coordinator for the black-oil simulator.
Definition: AdaptiveTimeStepping.hpp:93
double max_growth_
factor that limits the maximum growth of a time step
Definition: AdaptiveTimeStepping.hpp:435
double max_time_step_
maximal allowed time step size in days
Definition: AdaptiveTimeStepping.hpp:436
bool solver_verbose_
solver verbosity
Definition: AdaptiveTimeStepping.hpp:440
int solver_restart_max_
how many restart of solver are allowed
Definition: AdaptiveTimeStepping.hpp:439
double timestep_after_event_
suggested size of timestep after an event
Definition: AdaptiveTimeStepping.hpp:444
void init_(const UnitSystem &unitSystem)
Definition: AdaptiveTimeStepping_impl.hpp:429
void setSuggestedNextStep(const double x)
Set the suggested length for the next substep [s].
Definition: AdaptiveTimeStepping_impl.hpp:300
double suggestedNextStep() const
Definition: AdaptiveTimeStepping_impl.hpp:308
bool operator==(const AdaptiveTimeStepping< TypeTag > &rhs) const
Definition: AdaptiveTimeStepping_impl.hpp:140
static AdaptiveTimeStepping< TypeTag > serializationTestObjectSimple()
Definition: AdaptiveTimeStepping_impl.hpp:283
bool ignore_convergence_failure_
continue instead of stop when minimum time step is reached
Definition: AdaptiveTimeStepping.hpp:438
void serializeOp(Serializer &serializer)
Definition: AdaptiveTimeStepping_impl.hpp:212
void updateTUNING(double max_next_tstep, const Tuning &tuning)
Apply TUNING keyword parameters.
Definition: AdaptiveTimeStepping_impl.hpp:337
TimeStepControlType time_step_control_type_
type of time step control object
Definition: AdaptiveTimeStepping.hpp:431
const TimeStepControlInterface & timeStepControl() const
Definition: AdaptiveTimeStepping_impl.hpp:316
std::function< bool(double elapsed, double substep_length, int sub_step_number)> TuningUpdateCallback
Callback invoked at the start of each substep to apply TUNING, NEXTSTEP (via ACTIONX),...
Definition: AdaptiveTimeStepping.hpp:119
bool full_timestep_initially_
beginning with the size of the time step from data file
Definition: AdaptiveTimeStepping.hpp:443
SimulatorReport step(const SimulatorTimer &simulator_timer, Solver &solver, const bool is_event, const TuningUpdateCallback &tuning_updater)
Run one report step by orchestrating adaptive substepping.
Definition: AdaptiveTimeStepping_impl.hpp:197
double growth_factor_
factor to multiply time step when solver recovered from failed convergence
Definition: AdaptiveTimeStepping.hpp:434
double restart_factor_
factor to multiply time step with when solver fails to converge
Definition: AdaptiveTimeStepping.hpp:433
double min_time_step_
minimal allowed time step size before throwing
Definition: AdaptiveTimeStepping.hpp:437
void updateNEXTSTEP(double max_next_tstep)
Set suggested_next_timestep_ to max_next_tstep iff max_next_tstep > 0.
Definition: AdaptiveTimeStepping_impl.hpp:325
static AdaptiveTimeStepping< TypeTag > serializationTestObjectHardcoded()
Definition: AdaptiveTimeStepping_impl.hpp:259
TimeStepController time_step_control_
time step control object
Definition: AdaptiveTimeStepping.hpp:432
static AdaptiveTimeStepping< TypeTag > serializationTestObjectPIDIt()
Definition: AdaptiveTimeStepping_impl.hpp:275
double min_time_step_before_shutting_problematic_wells_
< shut problematic wells when time step size in days are less than this
Definition: AdaptiveTimeStepping.hpp:448
static AdaptiveTimeStepping< TypeTag > serializationTestObject3rdOrder()
Definition: AdaptiveTimeStepping_impl.hpp:291
SimulatorReport & report()
Definition: AdaptiveTimeStepping_impl.hpp:251
static AdaptiveTimeStepping< TypeTag > serializationTestObjectPID()
Definition: AdaptiveTimeStepping_impl.hpp:267
static void registerParameters()
Definition: AdaptiveTimeStepping_impl.hpp:186
bool use_newton_iteration_
use newton iteration count for adaptive time step control
Definition: AdaptiveTimeStepping.hpp:445
Definition: SimulatorTimer.hpp:38
Definition: TimeStepControlInterface.hpp:51
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition: parametersystem.hpp:191
void logTimer(const AdaptiveSimulatorTimer &substep_timer)
void registerAdaptiveParameters()
std::set< std::string > consistentlyFailingWells(const std::vector< StepReport > &sr, bool requireRepeatedFailures)
std::tuple< TimeStepControlType, std::unique_ptr< TimeStepControlInterface >, bool > createController(const UnitSystem &unitSystem)
Definition: blackoilbioeffectsmodules.hh:45
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
This file provides the infrastructure to retrieve run-time parameters.
static bool compare_gt_or_eq(double a, double b)
Determines if a is greater than b within the specified tolerance.
Definition: SimulatorReport.hpp:125