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
528// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()
529// It has to be called for each substep since TUNING might have been changed for next sub step due
530// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
531template<class TypeTag>
532template<class Solver>
533bool
534AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
535maybeUpdateTuning_(double elapsed, double substep_length, int sub_step_number) const
536{
537 return this->tuning_updater_(elapsed, substep_length, sub_step_number);
538}
539
540template<class TypeTag>
541template<class Solver>
542double
543AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
544maxTimeStep_() const
545{
546 return this->adaptive_time_stepping_.max_time_step_;
547}
548
549template <class TypeTag>
550template <class Solver>
551SimulatorReport
552AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
553runStepOriginal_()
554{
555 const auto elapsed = this->simulator_timer_.simulationTimeElapsed();
556 const auto original_time_step = this->simulator_timer_.currentStepLength();
557 const auto report_step = this->simulator_timer_.reportStepNum();
558 maybeUpdateTuning_(elapsed, suggestedNextTimestep_(), /*substep=*/0);
559 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(original_time_step);
560
561 AdaptiveSimulatorTimer substep_timer{
562 this->simulator_timer_.startDateTime(),
563 original_time_step,
564 elapsed,
565 suggestedNextTimestep_(),
566 report_step,
567 maxTimeStep_()
568 };
569 SubStepIteration<Solver> substepIteration{*this, substep_timer, original_time_step, /*final_step=*/true};
570 return substepIteration.run();
571}
572
573template <class TypeTag>
574template <class Solver>
575double
576AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
577suggestedNextTimestep_() const
578{
579 return this->adaptive_time_stepping_.suggestedNextStep();
580}
581
582
583#ifdef RESERVOIR_COUPLING_ENABLED
584// Pick the master's sync-step length for the next outer-loop iteration of
585// `runStepReservoirCouplingMaster_()`. Includes the chop against slave-
586// report dates and emits the user-visible log line. See the block comment
587// above `runStepReservoirCouplingMaster_()` for TSYNC/RSYNC mode definitions.
588//
589// TSYNC (default):
590// Each iteration is one master time step. Start from the master's
591// adaptive suggestion, capped at the max time step and at the time
592// remaining until the report-step end, then chop.
593//
594// RSYNC:
595// The outer loop iterates per chunk between slave-report boundaries.
596// Re-use the previous iteration's chopped value as the candidate for
597// this iteration (the caller seeds the first iteration with
598// `original_time_step`), then chop.
599template <class TypeTag>
600template <class Solver>
601double
602AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
603getRcMasterSyncStepLength_(double prev_step,
604 double current_time,
605 double step_end_time)
606{
607 const bool sync_at_report_steps = reservoirCouplingMaster_().syncAtReportSteps();
608 double current_step_length;
609 if (sync_at_report_steps) {
610 current_step_length = prev_step;
611 } else {
612 const double remaining = step_end_time - current_time;
613 current_step_length = std::min({suggestedNextTimestep_(), maxTimeStep_(), remaining});
614 // NOTE: The substep timer is later constructed (in the caller) with
615 // current_step_length as its span, and after construction,
616 // substep_timer.currentStepLength() should return the same as
617 // current_step_length. The timer's constructor calls
618 // provideTimeStepEstimate() with suggestedNextTimestep_() as the dt_
619 // estimate and current_step_length as the span. Since we took
620 // current_step_length = min(suggestedNextTimestep_(), maxTimeStep_(), remaining)
621 // here (and only shrink it further via maybeChopSubStep below), the
622 // estimate is always >= the span, so the timer's snap branch (see
623 // AdaptiveSimulatorTimer::provideTimeStepEstimate) fires and clamps
624 // dt_ to current_step_length.
625 }
626 current_step_length = reservoirCouplingMaster_().maybeChopSubStep(current_step_length, current_time);
627 auto num_active = reservoirCouplingMaster_().numActivatedSlaves();
628 OpmLog::info(fmt::format(
629 "\nChoosing next sync time{} between master and {} active slave {}: {:.2f} days",
630 sync_at_report_steps ? " (RSYNC)" : "",
631 num_active, (num_active == 1 ? "process" : "processes"),
632 current_step_length / unit::day
633 ));
634 return current_step_length;
635}
636
637template <class TypeTag>
638template <class Solver>
639ReservoirCouplingMaster<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
640AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
641reservoirCouplingMaster_()
642{
643 return *(this->solver_.model().simulator().reservoirCouplingMaster());
644}
645
646template <class TypeTag>
647template <class Solver>
648ReservoirCouplingSlave<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
649AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
650reservoirCouplingSlave_()
651{
652 return *(this->solver_.model().simulator().reservoirCouplingSlave());
653}
654
655// Description of the reservoir coupling master and slave substep loop
656// -------------------------------------------------------------------
657// The master and slave processes attempt to reach the end of the report step using a series of substeps
658// (also called timesteps). Each substep has an upper limit roughly determined by a combination of the
659// keywords TUNING (through TSINIT, TSMAXZ), NEXSTEP, WCYCLE, and the start of the next report step
660// (the last substep has to coincide with this time). Note that NEXTSTEP can be updated from an
661// ACTIONX keyword. Although this comment focuses on the maximum substep limit, there is also a lower
662// limit on the substep length, and the substep sizes are adjusted automatically (or retried) based on
663// the convergence behavior of the solver and other criteria.
664//
665// The master chooses a "sync step" — the span between successive master<->slave exchanges — via one
666// of two modes, selected by the `--rescoup-sync-at-report-steps` CLI flag:
667//
668// TSYNC (default, flag=false)
669// Each master time step is a sync step. Every outer iteration of the master substep
670// loop redistributes group-rate targets to the slaves. This gives fine sync granularity and
671// matches the reference simulator's coupling semantics.
672//
673// RSYNC (flag=true)
674// Each chunk between slave-report-step boundaries is one sync step. The master's own adaptive
675// timestepping runs inside the chunk, but slaves only hear back from the master once per chunk.
676// Provided as a developer escape hatch.
677//
678// In both modes the sync step is also limited so as not to overshoot the next slave report date
679// (`maybeChopSubStep`). A second, drift-based limit (GRUPMAST item 4 / RCMASTS) is parsed but not
680// yet honored — tracked as a follow-up PR.
681//
682// Per sync step:
683// - Master receives each activated slave's next report date (`receiveNextReportDateFromSlaves`).
684// - Master picks the sync-step length via `getRcMasterSyncStepLength_()` (mode-dependent) and
685// chops it so as not to overshoot any slave report date.
686// - Master sends the chosen sync-step length to each slave (`sendNextTimeStepToSlaves`).
687//
688// The slaves use the sync-step end as a fixed point — a "mini report step" — that they reach via one
689// or more of their own substeps. If the slave's current report step extends beyond the received
690// sync step, the slave loops waiting for the master to send further sync steps; the loop ends when
691// the received sync step coincides with the end of the slave's current report step.
692
693// Master-side rescoup outer loop. See the block comment above for TSYNC
694// and RSYNC mode definitions.
695//
696// TSYNC (default):
697// Each outer iteration is one master time step. The sync-step
698// length is chosen fresh each iteration by `getRcMasterSyncStepLength_()`
699// from the master's adaptive suggestion (capped at the max time step and
700// at the remaining report-step time), then chopped against slave-report
701// dates. `SubStepIteration::run()` is still used but normally handles
702// only solver-driven chops within a single master time step.
703//
704// RSYNC:
705// The outer loop iterates once per chunk between slave-report boundaries.
706// The sync-step length starts at the full report step and only shrinks
707// via `maybeChopSubStep`; the master's adaptive timestepping runs inside
708// `SubStepIteration` as solver-driven substeps within each chunk.
709template <class TypeTag>
710template <class Solver>
711SimulatorReport
712AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
713runStepReservoirCouplingMaster_()
714{
715 int iteration = 0;
716 const double original_time_step = this->simulator_timer_.currentStepLength();
717 double current_time{this->simulator_timer_.simulationTimeElapsed()};
718 double step_end_time = current_time + original_time_step;
719 // In RSYNC mode this variable persists across outer iterations and
720 // carries the previously-chopped sync span into the next iteration. In
721 // TSYNC mode it is overwritten each iteration by
722 // `getRcMasterSyncStepLength_()`; the initial value is only consumed by
723 // the first helper call in RSYNC mode (as the candidate span for that
724 // iteration).
725 auto current_step_length = original_time_step;
726 auto report_step_idx = this->simulator_timer_.currentStepNum();
727 if (report_step_idx == 0 && iteration == 0) {
728 reservoirCouplingMaster_().initTimeStepping();
729 }
730 SimulatorReport report;
731 // The master needs to know which slaves have activated before it can start the substep loop
732 reservoirCouplingMaster_().maybeReceiveActivationHandshakeFromSlaves(current_time);
733 while (true) {
734 reservoirCouplingMaster_().sendDontTerminateSignalToSlaves(); // Tell the slaves to keep running.
735 reservoirCouplingMaster_().receiveNextReportDateFromSlaves();
736 const bool start_of_report_step = (iteration == 0);
737 if (start_of_report_step) {
738 reservoirCouplingMaster_().initStartOfReportStep(report_step_idx);
739 maybeUpdateTuning_(current_time, suggestedNextTimestep_(), /*substep=*/0);
740 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(original_time_step);
741 }
742 current_step_length = getRcMasterSyncStepLength_(
743 current_step_length, current_time, step_end_time);
744 reservoirCouplingMaster_().sendNextTimeStepToSlaves(current_step_length);
745 AdaptiveSimulatorTimer substep_timer{
746 this->simulator_timer_.startDateTime(),
747 /*stepLength=*/current_step_length,
748 /*elapsedTime=*/current_time,
749 /*timeStepEstimate=*/suggestedNextTimestep_(),
750 /*reportStep=*/this->simulator_timer_.reportStepNum(),
751 maxTimeStep_()
752 };
753 const bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq(
754 current_time + current_step_length, step_end_time
755 );
756 // Mark this as the first substep of the "sync" timestep. This flag controls
757 // whether master-slave data exchange should occur in beginTimeStep() in the well model.
758 // It will be cleared after the first runSubStep_() call.
759 reservoirCouplingMaster_().setFirstSubstepOfSyncTimestep(true);
760 // After the first master substep completes, timeStepSucceeded() will
761 // block until slaves finish the sync step and send production data.
762 // This ensures correct summary output for all subsequent substeps.
763 reservoirCouplingMaster_().setNeedsSlaveDataReceive(true);
764 SubStepIteration<Solver> substepIteration{*this, substep_timer, current_step_length, final_step};
765 const auto sub_steps_report = substepIteration.run();
766 report += sub_steps_report;
767 current_time += current_step_length;
768 if (final_step) {
769 break;
770 }
771 iteration++;
772 }
773 return report;
774}
775
776template <class TypeTag>
777template <class Solver>
778SimulatorReport
779AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
780runStepReservoirCouplingSlave_()
781{
782 int iteration = 0;
783 const double original_time_step = this->simulator_timer_.currentStepLength();
784 double current_time{this->simulator_timer_.simulationTimeElapsed()};
785 double step_end_time = current_time + original_time_step;
786 SimulatorReport report;
787 auto report_step_idx = this->simulator_timer_.currentStepNum();
788 if (report_step_idx == 0 && iteration == 0) {
789 reservoirCouplingSlave_().initTimeStepping();
790 }
791 while (true) {
792 bool start_of_report_step = (iteration == 0);
793 if (reservoirCouplingSlave_().maybeReceiveTerminateSignalFromMaster()) {
794 // Call MPI_Comm_disconnect() to terminate the MPI communicator, etc..
795 break;
796 }
797 reservoirCouplingSlave_().sendNextReportDateToMasterProcess();
798 const auto timestep = reservoirCouplingSlave_().receiveNextTimeStepFromMaster();
799 if (start_of_report_step) {
800 maybeUpdateTuning_(current_time, suggestedNextTimestep_(), /*substep=*/0);
801 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(timestep);
802 }
803 AdaptiveSimulatorTimer substep_timer{
804 this->simulator_timer_.startDateTime(),
805 /*step_length=*/timestep,
806 /*elapsed_time=*/current_time,
807 /*time_step_estimate=*/suggestedNextTimestep_(),
808 this->simulator_timer_.reportStepNum(),
809 maxTimeStep_()
810 };
811 const bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq(
812 current_time + timestep, step_end_time
813 );
814 // Mark this as the first substep of the "sync" timestep. This flag controls
815 // whether master-slave data exchange should occur in beginTimeStep() in the well model.
816 // It will be cleared after the first runSubStep_() call.
817 reservoirCouplingSlave_().setFirstSubstepOfSyncTimestep(true);
818 SubStepIteration<Solver> substepIteration{*this, substep_timer, timestep, final_step};
819 const auto sub_steps_report = substepIteration.run();
820 report += sub_steps_report;
821 current_time += timestep;
822 if (final_step) {
823 break;
824 }
825 iteration++;
826 }
827 return report;
828}
829#endif // RESERVOIR_COUPLING_ENABLED
830
831
832/************************************************
833 * Private class SubStepIteration public methods
834 ************************************************/
835
836template<class TypeTag>
837template<class Solver>
838AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
839SubStepIteration(SubStepper<Solver>& substepper,
840 AdaptiveSimulatorTimer& substep_timer,
841 const double original_time_step,
842 bool final_step)
843 : substepper_{substepper}
844 , substep_timer_{substep_timer}
845 , original_time_step_{original_time_step}
846 , final_step_{final_step}
847 , adaptive_time_stepping_{substepper.getAdaptiveTimerStepper()}
848{
849}
850
851template <class TypeTag>
852template <class Solver>
853SimulatorReport
854AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
855run()
856{
857 auto& simulator = solver_().model().simulator();
858 auto& problem = simulator.problem();
859 // counter for solver restarts
860 int restarts = 0;
861 SimulatorReport report;
862
863 // sub step time loop
864 while (!this->substep_timer_.done()) {
865 // if we just chopped the timestep due to convergence i.e. restarts>0
866 // we dont want to update the next timestep based on Tuning
867 if (restarts == 0) {
868 maybeUpdateTuningAndTimeStep_();
869 }
870 const double dt = this->substep_timer_.currentStepLength();
871 if (timeStepVerbose_()) {
872 detail::logTimer(this->substep_timer_);
873 }
874
875 maybeUpdateLastSubstepOfSyncTimestep_(dt); // Needed for reservoir coupling
876 auto substep_report = runSubStep_();
877 markFirstSubStepAsFinished_(); // Needed for reservoir coupling
878
879 if (substep_report.converged || checkContinueOnUnconvergedSolution_(dt)) {
880 Dune::Timer perfTimer;
881 perfTimer.start();
882 // Pass substep to eclwriter for summary output
883 problem.setSubStepReport(substep_report);
884 auto& full_report = adaptive_time_stepping_.report();
885 full_report += substep_report;
886 problem.setSimulationReport(full_report);
887 problem.endTimeStep();
888 substep_report.pre_post_time += perfTimer.stop();
889
890 report += substep_report;
891
892 OPM_TIMEBLOCK(convergenceSucceeded);
893 ++this->substep_timer_; // advance by current dt
894
895 const int iterations = getNumIterations_(substep_report);
896 auto dt_estimate = timeStepControlComputeEstimate_(
897 dt, iterations, this->substep_timer_);
898
899 assert(dt_estimate > 0);
900 dt_estimate = maybeRestrictTimeStepGrowth_(dt, dt_estimate, restarts);
901 restarts = 0; // solver converged, reset restarts counter
902
903 maybeReportSubStep_(substep_report);
904 if (this->final_step_ && this->substep_timer_.done()) {
905 // if the time step is done we do not need to write it as this will be done
906 // by the simulator anyway.
907 }
908 else {
909 report.success.output_write_time += writeOutput_();
910 }
911
912 // set new time step length
913 setTimeStep_(dt_estimate);
914
915 report.success.converged = this->substep_timer_.done();
916 this->substep_timer_.setLastStepFailed(false);
917 }
918 else { // in case of no convergence or time step tolerance test failure
919 OPM_TIMEBLOCK(convergenceFailed);
920 report += substep_report;
921 this->substep_timer_.setLastStepFailed(true);
922 checkTimeStepMaxRestartLimit_(restarts);
923
924 double new_time_step = restartFactor_() * dt;
925 if (substep_report.time_step_rejected) {
926 const double tol = Parameters::Get<Parameters::TimeStepControlTolerance>();
927 const double safetyFactor = Parameters::Get<Parameters::TimeStepControlSafetyFactor>();
928 const double temp_time_step = std::sqrt(safetyFactor * tol / solver_().model().relativeChange()) * dt;
929 if (temp_time_step < dt) { // added in case suggested time step is not a reduction
930 new_time_step = temp_time_step;
931 }
932 }
933 checkTimeStepMinLimit_(new_time_step);
934 bool wells_shut = false;
935 if (new_time_step > minTimeStepBeforeClosingWells_()) {
936 chopTimeStep_(new_time_step);
937 } else {
938 wells_shut = chopTimeStepOrCloseFailingWells_(new_time_step);
939 }
940 if (wells_shut) {
941 setTimeStep_(dt); // retry the old timestep
942 }
943 else {
944 restarts++; // only increase if no wells were shut
945 }
946 }
947 problem.setNextTimeStepSize(this->substep_timer_.currentStepLength());
948 }
949 updateSuggestedNextStep_();
950 return report;
951}
952
953
954/************************************************
955 * Private class SubStepIteration private methods
956 ************************************************/
957
958
959template<class TypeTag>
960template<class Solver>
961bool
962AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
963checkContinueOnUnconvergedSolution_(double dt) const
964{
965 const bool continue_on_uncoverged_solution = ignoreConvergenceFailure_() && dt <= minTimeStep_();
966 if (continue_on_uncoverged_solution && solverVerbose_()) {
967 // NOTE: This method is only called if the solver failed to converge.
968 const auto msg = fmt::format(
969 "Solver failed to converge but timestep {} is smaller or equal to {}\n"
970 "which is the minimum threshold given by option --solver-min-time-step\n",
971 dt, minTimeStep_()
972 );
973 OpmLog::problem(msg);
974 }
975 return continue_on_uncoverged_solution;
976}
977
978template<class TypeTag>
979template<class Solver>
980void
981AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
982checkTimeStepMaxRestartLimit_(const int restarts) const
983{
984 // If we have restarted (i.e. cut the timestep) too
985 // many times, we have failed and throw an exception.
986 if (restarts >= solverRestartMax_()) {
987 const auto msg = fmt::format(
988 fmt::runtime("Solver failed to converge after cutting timestep {} times."), restarts
989 );
990 if (solverVerbose_()) {
991 OpmLog::error(msg);
992 }
993 // Use throw directly to prevent file and line
994 throw TimeSteppingBreakdown{msg};
995 }
996}
997
998template<class TypeTag>
999template<class Solver>
1000void
1001AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1002checkTimeStepMinLimit_(const double new_time_step) const
1003{
1004 using Meas = UnitSystem::measure;
1005 // If we have restarted (i.e. cut the timestep) too
1006 // much, we have failed and throw an exception.
1007 if (new_time_step < minTimeStep_()) {
1008 std::string msg = "Solver failed to converge after cutting timestep to ";
1009 if (Parameters::Get<Parameters::EnableTuning>()) {
1010 const UnitSystem& unit_system = solver_().model().simulator().vanguard().eclState().getDeckUnitSystem();
1011 msg += fmt::format(
1012 "{:.3E} {}\nwhich is the minimum threshold given by the TUNING keyword\n",
1013 unit_system.from_si(Meas::time, minTimeStep_()),
1014 unit_system.name(Meas::time)
1015 );
1016 }
1017 else {
1018 msg += fmt::format(
1019 "{:.3E} DAYS\nwhich is the minimum threshold given by option --solver-min-time-step\n",
1020 minTimeStep_() / 86400.0
1021 );
1022 }
1023 if (solverVerbose_()) {
1024 OpmLog::error(msg);
1025 }
1026 // Use throw directly to prevent file and line
1027 throw TimeSteppingBreakdown{msg};
1028 }
1029}
1030
1031template<class TypeTag>
1032template<class Solver>
1033void
1034AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1035chopTimeStep_(const double new_time_step)
1036{
1037 setTimeStep_(new_time_step);
1038 if (solverVerbose_()) {
1039 const auto msg = fmt::format(fmt::runtime("{}\nTimestep chopped to {} days\n"),
1040 this->cause_of_failure_,
1041 unit::convert::to(this->substep_timer_.currentStepLength(), unit::day));
1042 OpmLog::problem(msg);
1043 }
1044}
1045
1046template<class TypeTag>
1047template<class Solver>
1048bool
1049AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1050chopTimeStepOrCloseFailingWells_(const double new_time_step)
1051{
1052 bool wells_shut = false;
1053 // We are below the threshold, and will check if there are any
1054 // wells that fails repeatedly (that means that it fails in the last three steps)
1055 // we should close rather than chopping again.
1056 // If we already have chopped the timestep two times that is
1057 // new_time_step < minTimeStepBeforeClosingWells_()*restartFactor_()*restartFactor_()
1058 // We also shut wells that fails only on this step.
1059 const bool requireRepeatedFailures =
1060 new_time_step > (minTimeStepBeforeClosingWells_() * restartFactor_() * restartFactor_());
1061 const std::set<std::string> failing_wells =
1062 detail::consistentlyFailingWells(solver_().model().stepReports(), requireRepeatedFailures);
1063
1064 if (failing_wells.empty()) {
1065 // Found no wells to close, chop the timestep
1066 chopTimeStep_(new_time_step);
1067 } else {
1068 // Close all consistently failing wells that are not under group control
1069 std::vector<std::string> shut_wells;
1070 for (const auto& well : failing_wells) {
1071 const bool was_shut =
1072 solver_().model().wellModel().forceShutWellByName(well,
1073 this->substep_timer_.simulationTimeElapsed(),
1074 /*dont_shut_grup_wells =*/ true);
1075 if (was_shut) {
1076 shut_wells.push_back(well);
1077 }
1078 }
1079 // If no wells are closed we also try to shut wells under group control
1080 if (shut_wells.empty()) {
1081 for (const auto& well : failing_wells) {
1082 const bool was_shut =
1083 solver_().model().wellModel().forceShutWellByName(well,
1084 this->substep_timer_.simulationTimeElapsed(),
1085 /*dont_shut_grup_wells =*/ false);
1086 if (was_shut) {
1087 shut_wells.push_back(well);
1088 }
1089 }
1090 }
1091 // If still no wells are closed we must fall back to chopping again
1092 if (shut_wells.empty()) {
1093 chopTimeStep_(new_time_step);
1094 } else {
1095 wells_shut = true;
1096 if (solverVerbose_()) {
1097 const std::string msg =
1098 fmt::format(fmt::runtime("\nProblematic well(s) were shut: {}"
1099 "(retrying timestep)\n"),
1100 fmt::join(shut_wells, " "));
1101 OpmLog::problem(msg);
1102 }
1103 }
1104 }
1105 return wells_shut;
1106}
1107
1108template<class TypeTag>
1109template<class Solver>
1110boost::posix_time::ptime
1111AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1112currentDateTime_() const
1113{
1114 return simulatorTimer_().currentDateTime();
1115}
1116
1117template<class TypeTag>
1118template<class Solver>
1119int
1120AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1121getNumIterations_(const SimulatorReportSingle &substep_report) const
1122{
1123 if (useNewtonIteration_()) {
1124 return substep_report.total_newton_iterations;
1125 }
1126 else {
1127 return substep_report.total_linear_iterations;
1128 }
1129}
1130
1131template<class TypeTag>
1132template<class Solver>
1133double
1134AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1135growthFactor_() const
1136{
1137 return this->adaptive_time_stepping_.growth_factor_;
1138}
1139
1140template<class TypeTag>
1141template<class Solver>
1142bool
1143AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1144ignoreConvergenceFailure_() const
1145{
1146 return adaptive_time_stepping_.ignore_convergence_failure_;
1147}
1148
1149template<class TypeTag>
1150template<class Solver>
1151bool
1152AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1153isReservoirCouplingMaster_() const
1154{
1155 return this->substepper_.isReservoirCouplingMaster_();
1156}
1157
1158template<class TypeTag>
1159template<class Solver>
1160bool
1161AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1162isReservoirCouplingSlave_() const
1163{
1164 return this->substepper_.isReservoirCouplingSlave_();
1165}
1166
1167template<class TypeTag>
1168template<class Solver>
1169void
1170AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1171markFirstSubStepAsFinished_() const
1172{
1173#ifdef RESERVOIR_COUPLING_ENABLED
1174 // Clear the first-substep flag after the first runSubStep_() call.
1175 // This ensures that master-slave synchronization only happens once per sync timestep,
1176 // not on retry attempts after convergence-driven timestep chops.
1177 if (isReservoirCouplingMaster_()) {
1178 reservoirCouplingMaster_().setFirstSubstepOfSyncTimestep(false);
1179 }
1180 else if (isReservoirCouplingSlave_()) {
1181 reservoirCouplingSlave_().setFirstSubstepOfSyncTimestep(false);
1182 }
1183#endif
1184 return;
1185}
1186
1187template<class TypeTag>
1188template<class Solver>
1189double
1190AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1191maxGrowth_() const
1192{
1193 return this->adaptive_time_stepping_.max_growth_;
1194}
1195
1196template<class TypeTag>
1197template<class Solver>
1198void
1199AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1200maybeReportSubStep_(SimulatorReportSingle substep_report) const
1201{
1202 if (timeStepVerbose_()) {
1203 std::ostringstream ss;
1204 substep_report.reportStep(ss);
1205 OpmLog::info(ss.str());
1206 }
1207}
1208
1209template<class TypeTag>
1210template<class Solver>
1211double
1212AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1213maybeRestrictTimeStepGrowth_(const double dt, double dt_estimate, const int restarts) const
1214{
1215 // limit the growth of the timestep size by the growth factor
1216 dt_estimate = std::min(dt_estimate, double(maxGrowth_() * dt));
1217 assert(dt_estimate > 0);
1218 // further restrict time step size growth after convergence problems
1219 if (restarts > 0) {
1220 dt_estimate = std::min(growthFactor_() * dt, dt_estimate);
1221 }
1222
1223 return dt_estimate;
1224}
1225
1226
1227template<class TypeTag>
1228template<class Solver>
1229void
1230AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1231maybeUpdateLastSubstepOfSyncTimestep_([[maybe_unused]] const double dt)
1232{
1233#ifdef RESERVOIR_COUPLING_ENABLED
1234 // For reservoir coupling slaves: predict if this substep will complete
1235 // the sync timestep. If so, timeStepSucceeded() will send production
1236 // data to the master (which is blocking on receive after its first substep).
1237 // This is used for summary data synchronization between slaves and master.
1238 if (isReservoirCouplingSlave_()) {
1240 this->substep_timer_.simulationTimeElapsed() + dt,
1241 this->substep_timer_.totalTime()
1242 );
1243 reservoirCouplingSlave_().setLastSubstepOfSyncTimestep(is_last);
1244 }
1245#endif
1246}
1247
1248// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()
1249// It has to be called for each substep since TUNING might have been changed for next sub step due
1250// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
1251template<class TypeTag>
1252template<class Solver>
1253void
1254AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1255maybeUpdateTuningAndTimeStep_()
1256{
1257 // TODO: This function is currently only called if NEXTSTEP is activated from ACTIONX or
1258 // if the WCYCLE keyword needs to modify the current timestep. So this method should rather
1259 // be named maybeUpdateTimeStep_() or similar, since it should not update the tuning. However,
1260 // the current definition of the maybeUpdateTuning_() callback is actually calling
1261 // adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning) which is updating the tuning
1262 // see SimulatorFullyImplicitBlackoil::runStep() for more details.
1263 const auto old_value = suggestedNextTimestep_();
1264 if (this->substepper_.maybeUpdateTuning_(this->substep_timer_.simulationTimeElapsed(),
1265 this->substep_timer_.currentStepLength(),
1266 this->substep_timer_.currentStepNum()))
1267 {
1268 // Either NEXTSTEP and WCYCLE wants to change the current time step, but they cannot
1269 // change the current time step directly. Instead, they change the suggested next time step
1270 // by calling updateNEXTSTEP() via the maybeUpdateTuning() callback. We now need to update
1271 // the current time step to the new suggested time step and reset the suggested time step
1272 // to the old value.
1273 setTimeStep_(suggestedNextTimestep_());
1274 setSuggestedNextStep_(old_value);
1275 }
1276}
1277
1278template<class TypeTag>
1279template<class Solver>
1280double
1281AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1282minTimeStepBeforeClosingWells_() const
1283{
1284 return this->adaptive_time_stepping_.min_time_step_before_shutting_problematic_wells_;
1285}
1286
1287template<class TypeTag>
1288template<class Solver>
1289double
1290AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1291minTimeStep_() const
1292{
1293 return this->adaptive_time_stepping_.min_time_step_;
1294}
1295
1296#ifdef RESERVOIR_COUPLING_ENABLED
1297template<class TypeTag>
1298template<class Solver>
1299ReservoirCouplingMaster<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
1300AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1301reservoirCouplingMaster_() const
1302{
1303 return this->substepper_.reservoirCouplingMaster_();
1304}
1305
1306template<class TypeTag>
1307template<class Solver>
1308ReservoirCouplingSlave<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
1309AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1310reservoirCouplingSlave_() const
1311{
1312 return this->substepper_.reservoirCouplingSlave_();
1313}
1314#endif
1315
1316template<class TypeTag>
1317template<class Solver>
1318double
1319AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1320restartFactor_() const
1321{
1322 return this->adaptive_time_stepping_.restart_factor_;
1323}
1324
1325template<class TypeTag>
1326template<class Solver>
1327SimulatorReportSingle
1328AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1329runSubStep_()
1330{
1331 OPM_TIMEFUNCTION();
1332 SimulatorReportSingle substep_report;
1333
1334 auto handleFailure = [this, &substep_report]
1335 (const std::string& failure_reason, const std::exception& e, bool log_exception = true)
1336 {
1337 substep_report = solver_().failureReport();
1338 this->cause_of_failure_ = failure_reason;
1339 if (log_exception && solverVerbose_()) {
1340 OpmLog::debug(std::string("Caught Exception: ") + e.what());
1341 }
1342 };
1343
1344 try {
1345 substep_report = solver_().step(this->substep_timer_, &this->adaptive_time_stepping_.timeStepControl());
1346 if (solverVerbose_()) {
1347 // report number of linear iterations
1348 OpmLog::debug("Overall linear iterations used: "
1349 + std::to_string(substep_report.total_linear_iterations));
1350 }
1351 }
1352 catch (const TooManyIterations& e) {
1353 handleFailure("Solver convergence failure - Iteration limit reached", e);
1354 }
1355 catch (const TimeSteppingBreakdown& e) {
1356 handleFailure(e.what(), e);
1357 }
1358 catch (const ConvergenceMonitorFailure& e) {
1359 handleFailure("Convergence monitor failure", e, /*log_exception=*/false);
1360 }
1361 catch (const LinearSolverProblem& e) {
1362 handleFailure("Linear solver convergence failure", e);
1363 }
1364 catch (const NumericalProblem& e) {
1365 handleFailure("Solver convergence failure - Numerical problem encountered", e);
1366 }
1367 catch (const std::runtime_error& e) {
1368 handleFailure("Runtime error encountered", e);
1369 }
1370 catch (const Dune::ISTLError& e) {
1371 handleFailure("ISTL error - Time step too large", e);
1372 }
1373 catch (const Dune::MatrixBlockError& e) {
1374 handleFailure("Matrix block error", e);
1375 }
1376
1377 return substep_report;
1378}
1379
1380template<class TypeTag>
1381template<class Solver>
1382void
1383AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1384setTimeStep_(double dt_estimate)
1385{
1386 this->substep_timer_.provideTimeStepEstimate(dt_estimate);
1387}
1388
1389template<class TypeTag>
1390template<class Solver>
1391Solver&
1392AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1393solver_() const
1394{
1395 return this->substepper_.solver_;
1396}
1397
1398
1399template<class TypeTag>
1400template<class Solver>
1401int
1402AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1403solverRestartMax_() const
1404{
1405 return this->adaptive_time_stepping_.solver_restart_max_;
1406}
1407
1408template<class TypeTag>
1409template<class Solver>
1410void
1411AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1412setSuggestedNextStep_(double step)
1413{
1414 this->adaptive_time_stepping_.setSuggestedNextStep(step);
1415}
1416
1417template <class TypeTag>
1418template <class Solver>
1419const SimulatorTimer&
1420AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1421simulatorTimer_() const
1422{
1423 return this->substepper_.simulator_timer_;
1424}
1425
1426template <class TypeTag>
1427template <class Solver>
1428bool
1429AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1430solverVerbose_() const
1431{
1432 return this->adaptive_time_stepping_.solver_verbose_;
1433}
1434
1435template<class TypeTag>
1436template<class Solver>
1437boost::posix_time::ptime
1438AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1439startDateTime_() const
1440{
1441 return simulatorTimer_().startDateTime();
1442}
1443
1444template <class TypeTag>
1445template <class Solver>
1446double
1447AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1448suggestedNextTimestep_() const
1449{
1450 return this->adaptive_time_stepping_.suggestedNextStep();
1451}
1452
1453template <class TypeTag>
1454template <class Solver>
1455double
1456AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1457timeStepControlComputeEstimate_(const double dt, const int iterations,
1458 const AdaptiveSimulatorTimer& substepTimer) const
1459{
1460 // create object to compute the time error, simply forwards the call to the model
1461 const SolutionTimeErrorSolverWrapper<Solver> relative_change{solver_()};
1462 return this->adaptive_time_stepping_.time_step_control_->computeTimeStepSize(
1463 dt, iterations, relative_change, substepTimer);
1464}
1465
1466template <class TypeTag>
1467template <class Solver>
1468bool
1469AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1470timeStepVerbose_() const
1471{
1472 return this->adaptive_time_stepping_.timestep_verbose_;
1473}
1474
1475// The suggested time step is the stepsize that will be used as a first try for
1476// the next sub step. It is updated at the end of each substep. It can also be
1477// updated by the TUNING or NEXTSTEP keywords at the beginning of each report step or
1478// at the beginning of each substep by the ACTIONX keyword (via NEXTSTEP), this is
1479// done by the maybeUpdateTuning_() method which is called at the beginning of each substep
1480// (and the begginning of each report step). Note that the WCYCLE keyword can also update the
1481// suggested time step via the maybeUpdateTuning_() method.
1482template <class TypeTag>
1483template <class Solver>
1484void
1485AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1486updateSuggestedNextStep_()
1487{
1488 auto suggested_next_step = this->substep_timer_.currentStepLength();
1489 if (! std::isfinite(suggested_next_step)) { // check for NaN
1490 suggested_next_step = this->original_time_step_;
1491 }
1492 if (timeStepVerbose_()) {
1493 std::ostringstream ss;
1494 this->substep_timer_.report(ss);
1495 ss << "Suggested next step size = "
1496 << unit::convert::to(suggested_next_step, unit::day) << " (days)" << std::endl;
1497 OpmLog::debug(ss.str());
1498 }
1499 setSuggestedNextStep_(suggested_next_step);
1500}
1501
1502template <class TypeTag>
1503template <class Solver>
1504bool
1505AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1506useNewtonIteration_() const
1507{
1508 return this->adaptive_time_stepping_.use_newton_iteration_;
1509}
1510
1511template <class TypeTag>
1512template <class Solver>
1513double
1514AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1515writeOutput_() const
1516{
1517 time::StopWatch perf_timer;
1518 perf_timer.start();
1519 auto& problem = solver_().model().simulator().problem();
1520 problem.writeOutput(true);
1521 return perf_timer.secsSinceStart();
1522}
1523
1524/************************************************
1525 * Private class SolutionTimeErrorSolverWrapper
1526 * **********************************************/
1527
1528template<class TypeTag>
1529template<class Solver>
1530AdaptiveTimeStepping<TypeTag>::
1531SolutionTimeErrorSolverWrapper<Solver>::
1532SolutionTimeErrorSolverWrapper(const Solver& solver)
1533 : solver_{solver}
1534{}
1535
1536template<class TypeTag>
1537template<class Solver>
1538double AdaptiveTimeStepping<TypeTag>::SolutionTimeErrorSolverWrapper<Solver>::relativeChange() const
1539{
1540 // returns: || u^n+1 - u^n || / || u^n+1 ||
1541 return solver_.model().relativeChange();
1542}
1543
1544} // namespace Opm
1545
1546#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:427
double max_time_step_
maximal allowed time step size in days
Definition: AdaptiveTimeStepping.hpp:428
bool solver_verbose_
solver verbosity
Definition: AdaptiveTimeStepping.hpp:432
int solver_restart_max_
how many restart of solver are allowed
Definition: AdaptiveTimeStepping.hpp:431
double timestep_after_event_
suggested size of timestep after an event
Definition: AdaptiveTimeStepping.hpp:436
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:430
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:423
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:435
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:426
double restart_factor_
factor to multiply time step with when solver fails to converge
Definition: AdaptiveTimeStepping.hpp:425
double min_time_step_
minimal allowed time step size before throwing
Definition: AdaptiveTimeStepping.hpp:429
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:424
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:440
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:437
Definition: SimulatorTimer.hpp:39
Definition: TimeStepControlInterface.hpp:51
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition: parametersystem.hpp:187
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:122