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>
26#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
27#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
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
37#include <opm/grid/utility/StopWatch.hpp>
38
39#include <opm/input/eclipse/Schedule/Tuning.hpp>
40
41#include <opm/input/eclipse/Units/Units.hpp>
42#include <opm/input/eclipse/Units/UnitSystem.hpp>
43
44#include <opm/models/utils/parametersystem.hpp>
45
46#include <opm/simulators/timestepping/EclTimeSteppingParams.hpp>
47
48#include <algorithm>
49#include <cassert>
50#include <cmath>
51#include <sstream>
52#include <stdexcept>
53
54#include <boost/date_time/posix_time/posix_time.hpp>
55#include <fmt/format.h>
56#include <fmt/ranges.h>
57namespace Opm {
58/*********************************************
59 * Public methods of AdaptiveTimeStepping
60 * ******************************************/
61
62
64template<class TypeTag>
65AdaptiveTimeStepping<TypeTag>::
66AdaptiveTimeStepping(const UnitSystem& unit_system,
67 const SimulatorReport& report,
68 const double max_next_tstep,
69 const bool terminal_output
70)
71 : time_step_control_{}
72 , restart_factor_{Parameters::Get<Parameters::SolverRestartFactor<Scalar>>()} // 0.33
73 , growth_factor_{Parameters::Get<Parameters::SolverGrowthFactor<Scalar>>()} // 2.0
74 , max_growth_{Parameters::Get<Parameters::SolverMaxGrowth<Scalar>>()} // 3.0
75 , max_time_step_{
76 Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60} // 365.25
77 , min_time_step_{
78 unit_system.to_si(UnitSystem::measure::time,
79 Parameters::Get<Parameters::SolverMinTimeStep<Scalar>>())} // 1e-12;
80 , ignore_convergence_failure_{
81 Parameters::Get<Parameters::SolverContinueOnConvergenceFailure>()} // false;
82 , solver_restart_max_{Parameters::Get<Parameters::SolverMaxRestarts>()} // 10
83 , solver_verbose_{Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminal_output} // 2
84 , timestep_verbose_{Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminal_output} // 2
85 , suggested_next_timestep_{
86 (max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>()
87 : max_next_tstep) * 24 * 60 * 60} // 1.0
88 , full_timestep_initially_{Parameters::Get<Parameters::FullTimeStepInitially>()} // false
89 , timestep_after_event_{
90 Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60} // 1e30
91 , use_newton_iteration_{false}
92 , min_time_step_before_shutting_problematic_wells_{
93 Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day}
94 , report_(report)
95{
96 init_(unit_system);
97}
98
105template<class TypeTag>
106AdaptiveTimeStepping<TypeTag>::
107AdaptiveTimeStepping(double max_next_tstep,
108 const Tuning& tuning,
109 const UnitSystem& unit_system,
110 const SimulatorReport& report,
111 const bool terminal_output
112)
113 : time_step_control_{}
114 , restart_factor_{tuning.TSFCNV}
115 , growth_factor_{tuning.TFDIFF}
116 , max_growth_{tuning.TSFMAX}
117 , max_time_step_{tuning.TSMAXZ} // 365.25
118 , min_time_step_{tuning.TSMINZ} // 0.1;
119 , ignore_convergence_failure_{true}
120 , solver_restart_max_{Parameters::Get<Parameters::SolverMaxRestarts>()} // 10
121 , solver_verbose_{Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminal_output} // 2
122 , timestep_verbose_{Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminal_output} // 2
123 , suggested_next_timestep_{
124 max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() * 24 * 60 * 60
125 : max_next_tstep} // 1.0
126 , full_timestep_initially_{Parameters::Get<Parameters::FullTimeStepInitially>()} // false
127 , timestep_after_event_{tuning.TMAXWC} // 1e30
128 , use_newton_iteration_{false}
129 , min_time_step_before_shutting_problematic_wells_{
130 Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day}
131 , report_(report)
132{
133 init_(unit_system);
134}
135
136template<class TypeTag>
137bool
138AdaptiveTimeStepping<TypeTag>::
139operator==(const AdaptiveTimeStepping<TypeTag>& rhs) const
140{
141 if (this->time_step_control_type_ != rhs.time_step_control_type_ ||
142 (this->time_step_control_ && !rhs.time_step_control_) ||
143 (!this->time_step_control_ && rhs.time_step_control_)) {
144 return false;
145 }
146
147 bool result = false;
148 switch (this->time_step_control_type_) {
149 case TimeStepControlType::HardCodedTimeStep:
150 result = castAndComp<HardcodedTimeStepControl>(rhs);
151 break;
152 case TimeStepControlType::PIDAndIterationCount:
153 result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
154 break;
155 case TimeStepControlType::SimpleIterationCount:
156 result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
157 break;
158 case TimeStepControlType::PID:
159 result = castAndComp<PIDTimeStepControl>(rhs);
160 break;
161 case TimeStepControlType::General3rdOrder:
162 result = castAndComp<General3rdOrderController>(rhs);
163 break;
164 }
165
166 return result &&
167 this->restart_factor_ == rhs.restart_factor_ &&
168 this->growth_factor_ == rhs.growth_factor_ &&
169 this->max_growth_ == rhs.max_growth_ &&
170 this->max_time_step_ == rhs.max_time_step_ &&
171 this->min_time_step_ == rhs.min_time_step_ &&
172 this->ignore_convergence_failure_ == rhs.ignore_convergence_failure_ &&
173 this->solver_restart_max_== rhs.solver_restart_max_ &&
174 this->solver_verbose_ == rhs.solver_verbose_ &&
175 this->full_timestep_initially_ == rhs.full_timestep_initially_ &&
176 this->timestep_after_event_ == rhs.timestep_after_event_ &&
177 this->use_newton_iteration_ == rhs.use_newton_iteration_ &&
178 this->min_time_step_before_shutting_problematic_wells_ ==
179 rhs.min_time_step_before_shutting_problematic_wells_;
180}
181
182template<class TypeTag>
183void
184AdaptiveTimeStepping<TypeTag>::
185registerParameters()
186{
187 registerEclTimeSteppingParameters<Scalar>();
188 detail::registerAdaptiveParameters();
189}
190
191// See Doxygen comment on the declaration in AdaptiveTimeStepping.hpp.
192template<class TypeTag>
193template <class Solver>
194SimulatorReport
195AdaptiveTimeStepping<TypeTag>::
196step(const SimulatorTimer& simulator_timer,
197 Solver& solver,
198 const bool is_event,
199 const TuningUpdateCallback& tuning_updater)
200{
201 SubStepper<Solver> sub_stepper{
202 *this, simulator_timer, solver, is_event, tuning_updater,
203 };
204 return sub_stepper.run();
205}
206
207template<class TypeTag>
208template<class Serializer>
209void
210AdaptiveTimeStepping<TypeTag>::
211serializeOp(Serializer& serializer)
212{
213 serializer(this->time_step_control_type_);
214 switch (this->time_step_control_type_) {
215 case TimeStepControlType::HardCodedTimeStep:
216 allocAndSerialize<HardcodedTimeStepControl>(serializer);
217 break;
218 case TimeStepControlType::PIDAndIterationCount:
219 allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
220 break;
221 case TimeStepControlType::SimpleIterationCount:
222 allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
223 break;
224 case TimeStepControlType::PID:
225 allocAndSerialize<PIDTimeStepControl>(serializer);
226 break;
227 case TimeStepControlType::General3rdOrder:
228 allocAndSerialize<General3rdOrderController>(serializer);
229 break;
230 }
231 serializer(this->restart_factor_);
232 serializer(this->growth_factor_);
233 serializer(this->max_growth_);
234 serializer(this->max_time_step_);
235 serializer(this->min_time_step_);
236 serializer(this->ignore_convergence_failure_);
237 serializer(this->solver_restart_max_);
238 serializer(this->solver_verbose_);
239 serializer(this->timestep_verbose_);
240 serializer(this->suggested_next_timestep_);
241 serializer(this->full_timestep_initially_);
242 serializer(this->timestep_after_event_);
243 serializer(this->use_newton_iteration_);
244 serializer(this->min_time_step_before_shutting_problematic_wells_);
245}
246
247template<class TypeTag>
248SimulatorReport&
249AdaptiveTimeStepping<TypeTag>::
250report()
251{
252 return report_;
253}
254
255template<class TypeTag>
256AdaptiveTimeStepping<TypeTag>
257AdaptiveTimeStepping<TypeTag>::
258serializationTestObjectHardcoded()
259{
260 return serializationTestObject_<HardcodedTimeStepControl>();
261}
262
263template<class TypeTag>
264AdaptiveTimeStepping<TypeTag>
265AdaptiveTimeStepping<TypeTag>::
266serializationTestObjectPID()
267{
268 return serializationTestObject_<PIDTimeStepControl>();
269}
270
271template<class TypeTag>
272AdaptiveTimeStepping<TypeTag>
273AdaptiveTimeStepping<TypeTag>::
274serializationTestObjectPIDIt()
275{
276 return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
277}
278
279template<class TypeTag>
280AdaptiveTimeStepping<TypeTag>
281AdaptiveTimeStepping<TypeTag>::
282serializationTestObjectSimple()
283{
284 return serializationTestObject_<SimpleIterationCountTimeStepControl>();
285}
286
287template<class TypeTag>
288AdaptiveTimeStepping<TypeTag>
289AdaptiveTimeStepping<TypeTag>::
290serializationTestObject3rdOrder()
291{
292 return serializationTestObject_<General3rdOrderController>();
293}
294
295
296template<class TypeTag>
297void
298AdaptiveTimeStepping<TypeTag>::
299setSuggestedNextStep(const double x)
300{
301 this->suggested_next_timestep_ = x;
302}
303
304template<class TypeTag>
305double
306AdaptiveTimeStepping<TypeTag>::
307suggestedNextStep() const
308{
309 return this->suggested_next_timestep_;
310}
311
312template<class TypeTag>
313const TimeStepControlInterface&
314AdaptiveTimeStepping<TypeTag>::
315timeStepControl() const
316{
317 return *this->time_step_control_;
318}
319
320
321template<class TypeTag>
322void
323AdaptiveTimeStepping<TypeTag>::
324updateNEXTSTEP(double max_next_tstep)
325{
326 // \Note Only update next suggested step if TSINIT was explicitly
327 // set in TUNING or NEXTSTEP is active.
328 if (max_next_tstep > 0) {
329 this->suggested_next_timestep_ = max_next_tstep;
330 }
331}
332
333template<class TypeTag>
334void
335AdaptiveTimeStepping<TypeTag>::
336updateTUNING(double max_next_tstep, const Tuning& tuning)
337{
338 this->restart_factor_ = tuning.TSFCNV;
339 this->growth_factor_ = tuning.TFDIFF;
340 this->max_growth_ = tuning.TSFMAX;
341 this->max_time_step_ = tuning.TSMAXZ;
342 updateNEXTSTEP(max_next_tstep);
343 this->timestep_after_event_ = tuning.TMAXWC;
344}
345
346/*********************************************
347 * Private methods of AdaptiveTimeStepping
348 * ******************************************/
349
350template<class TypeTag>
351template<class T, class Serializer>
352void
353AdaptiveTimeStepping<TypeTag>::
354allocAndSerialize(Serializer& serializer)
355{
356 if (!serializer.isSerializing()) {
357 this->time_step_control_ = std::make_unique<T>();
358 }
359 serializer(*static_cast<T*>(this->time_step_control_.get()));
360}
361
362template<class TypeTag>
363template<class T>
364bool
365AdaptiveTimeStepping<TypeTag>::
366castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const
367{
368 const T* lhs = static_cast<const T*>(this->time_step_control_.get());
369 const T* rhs = static_cast<const T*>(Rhs.time_step_control_.get());
370 return *lhs == *rhs;
371}
372
373template<class TypeTag>
374void
375AdaptiveTimeStepping<TypeTag>::
376maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step,
377 bool is_event)
378{
379 // init last time step as a fraction of the given time step
380 if (this->suggested_next_timestep_ < 0) {
381 this->suggested_next_timestep_ = this->restart_factor_ * original_time_step;
382 }
383
384 if (this->full_timestep_initially_) {
385 this->suggested_next_timestep_ = original_time_step;
386 }
387
388 // use seperate time step after event
389 if (is_event && this->timestep_after_event_ > 0) {
390 this->suggested_next_timestep_ = this->timestep_after_event_;
391 }
392}
393
394template<class TypeTag>
395template<class Controller>
396AdaptiveTimeStepping<TypeTag>
397AdaptiveTimeStepping<TypeTag>::
398serializationTestObject_()
399{
400 AdaptiveTimeStepping<TypeTag> result;
401
402 result.restart_factor_ = 1.0;
403 result.growth_factor_ = 2.0;
404 result.max_growth_ = 3.0;
405 result.max_time_step_ = 4.0;
406 result.min_time_step_ = 5.0;
407 result.ignore_convergence_failure_ = true;
408 result.solver_restart_max_ = 6;
409 result.solver_verbose_ = true;
410 result.timestep_verbose_ = true;
411 result.suggested_next_timestep_ = 7.0;
412 result.full_timestep_initially_ = true;
413 result.use_newton_iteration_ = true;
414 result.min_time_step_before_shutting_problematic_wells_ = 9.0;
415 result.time_step_control_type_ = Controller::Type;
416 result.time_step_control_ =
417 std::make_unique<Controller>(Controller::serializationTestObject());
418
419 return result;
420}
421
422/*********************************************
423 * Protected methods of AdaptiveTimeStepping
424 * ******************************************/
425
426template<class TypeTag>
427void AdaptiveTimeStepping<TypeTag>::
428init_(const UnitSystem& unitSystem)
429{
430 std::tie(time_step_control_type_,
431 time_step_control_,
432 use_newton_iteration_) = detail::createController(unitSystem);
433 // make sure growth factor is something reasonable
434 if (this->growth_factor_ < 1.0) {
435 OPM_THROW(std::runtime_error,
436 "Growth factor cannot be less than 1.");
437 }
438}
439
440
441
442/************************************************
443 * Private class SubStepper public methods
444 ************************************************/
445
446template<class TypeTag>
447template<class Solver>
448AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
449SubStepper(AdaptiveTimeStepping<TypeTag>& adaptive_time_stepping,
450 const SimulatorTimer& simulator_timer,
451 Solver& solver,
452 const bool is_event,
453 const TuningUpdateCallback& tuning_updater)
454 : adaptive_time_stepping_{adaptive_time_stepping}
455 , simulator_timer_{simulator_timer}
456 , solver_{solver}
457 , is_event_{is_event}
458 , tuning_updater_{tuning_updater}
459{
460}
461
462template<class TypeTag>
463template<class Solver>
464AdaptiveTimeStepping<TypeTag>&
465AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
466getAdaptiveTimerStepper()
467{
468 return adaptive_time_stepping_;
469}
470
471template<class TypeTag>
472template<class Solver>
473SimulatorReport
474AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
475run()
476{
477#ifdef RESERVOIR_COUPLING_ENABLED
478 if (isReservoirCouplingSlave_() && reservoirCouplingSlave_().activated()) {
479 return runStepReservoirCouplingSlave_();
480 }
481 else if (isReservoirCouplingMaster_() && reservoirCouplingMaster_().activated()) {
482 return runStepReservoirCouplingMaster_();
483 }
484 else {
485 return runStepOriginal_();
486 }
487#else
488 return runStepOriginal_();
489#endif
490}
491
492/************************************************
493 * Private class SubStepper private methods
494 ************************************************/
495
496#ifdef RESERVOIR_COUPLING_ENABLED
497template<class TypeTag>
498template<class Solver>
499bool
500AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
501isReservoirCouplingMaster_() const
502{
503 return this->solver_.model().simulator().reservoirCouplingMaster() != nullptr;
504}
505
506template<class TypeTag>
507template<class Solver>
508bool
509AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
510isReservoirCouplingSlave_() const
511{
512 return this->solver_.model().simulator().reservoirCouplingSlave() != nullptr;
513}
514#endif
515
516template<class TypeTag>
517template<class Solver>
518void
519AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
520maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step)
521{
522 this->adaptive_time_stepping_.maybeModifySuggestedTimeStepAtBeginningOfReportStep_(
523 original_time_step, this->is_event_
524 );
525}
526
527// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()
528// It has to be called for each substep since TUNING might have been changed for next sub step due
529// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
530template<class TypeTag>
531template<class Solver>
532bool
533AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
534maybeUpdateTuning_(double elapsed, double substep_length, int sub_step_number) const
535{
536 return this->tuning_updater_(elapsed, substep_length, sub_step_number);
537}
538
539template<class TypeTag>
540template<class Solver>
541double
542AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
543maxTimeStep_() const
544{
545 return this->adaptive_time_stepping_.max_time_step_;
546}
547
548template <class TypeTag>
549template <class Solver>
550SimulatorReport
551AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
552runStepOriginal_()
553{
554 const auto elapsed = this->simulator_timer_.simulationTimeElapsed();
555 const auto original_time_step = this->simulator_timer_.currentStepLength();
556 const auto report_step = this->simulator_timer_.reportStepNum();
557 maybeUpdateTuning_(elapsed, suggestedNextTimestep_(), /*substep=*/0);
558 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(original_time_step);
559
560 AdaptiveSimulatorTimer substep_timer{
561 this->simulator_timer_.startDateTime(),
562 original_time_step,
563 elapsed,
564 suggestedNextTimestep_(),
565 report_step,
566 maxTimeStep_()
567 };
568 SubStepIteration<Solver> substepIteration{*this, substep_timer, original_time_step, /*final_step=*/true};
569 return substepIteration.run();
570}
571
572#ifdef RESERVOIR_COUPLING_ENABLED
573template <class TypeTag>
574template <class Solver>
575ReservoirCouplingMaster<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
576AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
577reservoirCouplingMaster_()
578{
579 return *(this->solver_.model().simulator().reservoirCouplingMaster());
580}
581#endif
582
583#ifdef RESERVOIR_COUPLING_ENABLED
584template <class TypeTag>
585template <class Solver>
586ReservoirCouplingSlave<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
587AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
588reservoirCouplingSlave_()
589{
590 return *(this->solver_.model().simulator().reservoirCouplingSlave());
591}
592#endif
593
594#ifdef RESERVOIR_COUPLING_ENABLED
595// Description of the reservoir coupling master and slave substep loop
596// -------------------------------------------------------------------
597// The master and slave processes attempts to reach the end of the report step using a series of substeps
598// (also called timesteps). Each substep have an upper limit that is roughly determined by a combination
599// of the keywords TUNING (through the TSINIT, TSMAXZ values), NEXSTEP, WCYCLE, and the start of the
600// next report step (the last substep needs to coincide with this time). Note that NEXTSTEP can be
601// updated from an ACTIONX keyword. Although, this comment focuses on the maximum substep limit,
602// note that there is also a lower limit on the substep length. And the substep sizes will be adjusted
603// automatically (or retried) based on the convergence behavior of the solver and other criteria.
604//
605// When using reservoir coupling, the upper limit on each substep is further limited to: a) not overshoot
606// next report date of a slave reservoir, and b) to keep the flow rate of the slave groups within
607// certain limits. To determine this potential further limitation on the substep, the following procedure
608// is used at the beginning of each master substep:
609// - First, the slaves sends the master the date of their next report step
610// - The master receives the date of the next report step from the slaves
611// - If necessary, the master computes a modified substep length based on the received dates
612// TODO: explain how the substep is limited to keep the flow rate of the slave groups within certain
613// limits. Since this is not implemented yet, this part of the procedure is not explained here.
614//
615// Then, after the master has determined the substep length using the above procedure, it sends it
616// to the slaves. The slaves will now use the end of this substep as a fixed point (like a mini report
617// step), that they will try to reach using a single substep or multiple substeps. The end of the
618// substep received from the master will either conincide with the end of its current report step, or
619// it will end before (it cannot not end after since the master has already adjusted the step to not
620// exceed any report time in a slave). If the end of this substep is earlier in time than its next report
621// date, the slave will enter a loop and wait for the master to send it a new substep.
622// The loop is finished when the end of the received time step conincides with the end of its current
623// report step.
624
625template <class TypeTag>
626template <class Solver>
627SimulatorReport
628AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
629runStepReservoirCouplingMaster_()
630{
631 int iteration = 0;
632 const double original_time_step = this->simulator_timer_.currentStepLength();
633 double current_time{this->simulator_timer_.simulationTimeElapsed()};
634 double step_end_time = current_time + original_time_step;
635 auto current_step_length = original_time_step;
636 auto report_step_idx = this->simulator_timer_.currentStepNum();
637 if (report_step_idx == 0 && iteration == 0) {
638 reservoirCouplingMaster_().initTimeStepping();
639 }
640 SimulatorReport report;
641 // The master needs to know which slaves have activated before it can start the substep loop
642 reservoirCouplingMaster_().maybeReceiveActivationHandshakeFromSlaves(current_time);
643 while (true) {
644 reservoirCouplingMaster_().sendDontTerminateSignalToSlaves(); // Tell the slaves to keep running.
645 reservoirCouplingMaster_().receiveNextReportDateFromSlaves();
646 bool start_of_report_step = (iteration == 0);
647 if (start_of_report_step) {
648 reservoirCouplingMaster_().initStartOfReportStep(report_step_idx);
649 }
650 current_step_length = reservoirCouplingMaster_().maybeChopSubStep(
651 current_step_length, current_time);
652 auto num_active = reservoirCouplingMaster_().numActivatedSlaves();
653 OpmLog::info(fmt::format(
654 "\nChoosing next sync time between master and {} active slave {}: {:.2f} days",
655 num_active, (num_active == 1 ? "process" : "processes"),
656 current_step_length / unit::day
657 ));
658 reservoirCouplingMaster_().sendNextTimeStepToSlaves(current_step_length);
659 if (start_of_report_step) {
660 maybeUpdateTuning_(current_time, suggestedNextTimestep_(), /*substep=*/0);
661 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(current_step_length);
662 }
663 AdaptiveSimulatorTimer substep_timer{
664 this->simulator_timer_.startDateTime(),
665 /*stepLength=*/current_step_length,
666 /*elapsedTime=*/current_time,
667 /*timeStepEstimate=*/suggestedNextTimestep_(),
668 /*reportStep=*/this->simulator_timer_.reportStepNum(),
669 maxTimeStep_()
670 };
671 const bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq(
672 current_time + current_step_length, step_end_time
673 );
674 // Mark this as the first substep of the "sync" timestep. This flag controls
675 // whether master-slave data exchange should occur in beginTimeStep() in the well model.
676 // It will be cleared after the first runSubStep_() call.
677 reservoirCouplingMaster_().setFirstSubstepOfSyncTimestep(true);
678 // After the first master substep completes, timeStepSucceeded() will
679 // block until slaves finish the sync step and send production data.
680 // This ensures correct summary output for all subsequent substeps.
681 reservoirCouplingMaster_().setNeedsSlaveDataReceive(true);
682 SubStepIteration<Solver> substepIteration{*this, substep_timer, current_step_length, final_step};
683 const auto sub_steps_report = substepIteration.run();
684 report += sub_steps_report;
685 current_time += current_step_length;
686 if (final_step) {
687 break;
688 }
689 iteration++;
690 }
691 return report;
692}
693#endif
694
695#ifdef RESERVOIR_COUPLING_ENABLED
696template <class TypeTag>
697template <class Solver>
698SimulatorReport
699AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
700runStepReservoirCouplingSlave_()
701{
702 int iteration = 0;
703 const double original_time_step = this->simulator_timer_.currentStepLength();
704 double current_time{this->simulator_timer_.simulationTimeElapsed()};
705 double step_end_time = current_time + original_time_step;
706 SimulatorReport report;
707 auto report_step_idx = this->simulator_timer_.currentStepNum();
708 if (report_step_idx == 0 && iteration == 0) {
709 reservoirCouplingSlave_().initTimeStepping();
710 }
711 while (true) {
712 bool start_of_report_step = (iteration == 0);
713 if (reservoirCouplingSlave_().maybeReceiveTerminateSignalFromMaster()) {
714 // Call MPI_Comm_disconnect() to terminate the MPI communicator, etc..
715 break;
716 }
717 reservoirCouplingSlave_().sendNextReportDateToMasterProcess();
718 const auto timestep = reservoirCouplingSlave_().receiveNextTimeStepFromMaster();
719 if (start_of_report_step) {
720 maybeUpdateTuning_(current_time, suggestedNextTimestep_(), /*substep=*/0);
721 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(timestep);
722 }
723 AdaptiveSimulatorTimer substep_timer{
724 this->simulator_timer_.startDateTime(),
725 /*step_length=*/timestep,
726 /*elapsed_time=*/current_time,
727 /*time_step_estimate=*/suggestedNextTimestep_(),
728 this->simulator_timer_.reportStepNum(),
729 maxTimeStep_()
730 };
731 const bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq(
732 current_time + timestep, step_end_time
733 );
734 // Mark this as the first substep of the "sync" timestep. This flag controls
735 // whether master-slave data exchange should occur in beginTimeStep() in the well model.
736 // It will be cleared after the first runSubStep_() call.
737 reservoirCouplingSlave_().setFirstSubstepOfSyncTimestep(true);
738 SubStepIteration<Solver> substepIteration{*this, substep_timer, timestep, final_step};
739 const auto sub_steps_report = substepIteration.run();
740 report += sub_steps_report;
741 current_time += timestep;
742 if (final_step) {
743 break;
744 }
745 iteration++;
746 }
747 return report;
748}
749
750#endif
751
752template <class TypeTag>
753template <class Solver>
754double
755AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
756suggestedNextTimestep_() const
757{
758 return this->adaptive_time_stepping_.suggestedNextStep();
759}
760
761
762
763/************************************************
764 * Private class SubStepIteration public methods
765 ************************************************/
766
767template<class TypeTag>
768template<class Solver>
769AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
770SubStepIteration(SubStepper<Solver>& substepper,
771 AdaptiveSimulatorTimer& substep_timer,
772 const double original_time_step,
773 bool final_step)
774 : substepper_{substepper}
775 , substep_timer_{substep_timer}
776 , original_time_step_{original_time_step}
777 , final_step_{final_step}
778 , adaptive_time_stepping_{substepper.getAdaptiveTimerStepper()}
779{
780}
781
782template <class TypeTag>
783template <class Solver>
784SimulatorReport
785AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
786run()
787{
788 auto& simulator = solver_().model().simulator();
789 auto& problem = simulator.problem();
790 // counter for solver restarts
791 int restarts = 0;
792 SimulatorReport report;
793
794 // sub step time loop
795 while (!this->substep_timer_.done()) {
796 // if we just chopped the timestep due to convergence i.e. restarts>0
797 // we dont want to update the next timestep based on Tuning
798 if (restarts == 0) {
799 maybeUpdateTuningAndTimeStep_();
800 }
801 const double dt = this->substep_timer_.currentStepLength();
802 if (timeStepVerbose_()) {
803 detail::logTimer(this->substep_timer_);
804 }
805
806 maybeUpdateLastSubstepOfSyncTimestep_(dt); // Needed for reservoir coupling
807 auto substep_report = runSubStep_();
808 markFirstSubStepAsFinished_(); // Needed for reservoir coupling
809
810 if (substep_report.converged || checkContinueOnUnconvergedSolution_(dt)) {
811 Dune::Timer perfTimer;
812 perfTimer.start();
813 // Pass substep to eclwriter for summary output
814 problem.setSubStepReport(substep_report);
815 auto& full_report = adaptive_time_stepping_.report();
816 full_report += substep_report;
817 problem.setSimulationReport(full_report);
818 problem.endTimeStep();
819 substep_report.pre_post_time += perfTimer.stop();
820
821 report += substep_report;
822
823 ++this->substep_timer_; // advance by current dt
824
825 const int iterations = getNumIterations_(substep_report);
826 auto dt_estimate = timeStepControlComputeEstimate_(
827 dt, iterations, this->substep_timer_);
828
829 assert(dt_estimate > 0);
830 dt_estimate = maybeRestrictTimeStepGrowth_(dt, dt_estimate, restarts);
831 restarts = 0; // solver converged, reset restarts counter
832
833 maybeReportSubStep_(substep_report);
834 if (this->final_step_ && this->substep_timer_.done()) {
835 // if the time step is done we do not need to write it as this will be done
836 // by the simulator anyway.
837 }
838 else {
839 report.success.output_write_time += writeOutput_();
840 }
841
842 // set new time step length
843 setTimeStep_(dt_estimate);
844
845 report.success.converged = this->substep_timer_.done();
846 this->substep_timer_.setLastStepFailed(false);
847 }
848 else { // in case of no convergence or time step tolerance test failure
849 report += substep_report;
850 this->substep_timer_.setLastStepFailed(true);
851 checkTimeStepMaxRestartLimit_(restarts);
852
853 double new_time_step = restartFactor_() * dt;
854 if (substep_report.time_step_rejected) {
855 const double tol = Parameters::Get<Parameters::TimeStepControlTolerance>();
856 const double safetyFactor = Parameters::Get<Parameters::TimeStepControlSafetyFactor>();
857 const double temp_time_step = std::sqrt(safetyFactor * tol / solver_().model().relativeChange()) * dt;
858 if (temp_time_step < dt) { // added in case suggested time step is not a reduction
859 new_time_step = temp_time_step;
860 }
861 }
862 checkTimeStepMinLimit_(new_time_step);
863 bool wells_shut = false;
864 if (new_time_step > minTimeStepBeforeClosingWells_()) {
865 chopTimeStep_(new_time_step);
866 } else {
867 wells_shut = chopTimeStepOrCloseFailingWells_(new_time_step);
868 }
869 if (wells_shut) {
870 setTimeStep_(dt); // retry the old timestep
871 }
872 else {
873 restarts++; // only increase if no wells were shut
874 }
875 }
876 problem.setNextTimeStepSize(this->substep_timer_.currentStepLength());
877 }
878 updateSuggestedNextStep_();
879 return report;
880}
881
882
883/************************************************
884 * Private class SubStepIteration private methods
885 ************************************************/
886
887
888template<class TypeTag>
889template<class Solver>
890bool
891AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
892checkContinueOnUnconvergedSolution_(double dt) const
893{
894 const bool continue_on_uncoverged_solution = ignoreConvergenceFailure_() && dt <= minTimeStep_();
895 if (continue_on_uncoverged_solution && solverVerbose_()) {
896 // NOTE: This method is only called if the solver failed to converge.
897 const auto msg = fmt::format(
898 "Solver failed to converge but timestep {} is smaller or equal to {}\n"
899 "which is the minimum threshold given by option --solver-min-time-step\n",
900 dt, minTimeStep_()
901 );
902 OpmLog::problem(msg);
903 }
904 return continue_on_uncoverged_solution;
905}
906
907template<class TypeTag>
908template<class Solver>
909void
910AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
911checkTimeStepMaxRestartLimit_(const int restarts) const
912{
913 // If we have restarted (i.e. cut the timestep) too
914 // many times, we have failed and throw an exception.
915 if (restarts >= solverRestartMax_()) {
916 const auto msg = fmt::format(
917 fmt::runtime("Solver failed to converge after cutting timestep {} times."), restarts
918 );
919 if (solverVerbose_()) {
920 OpmLog::error(msg);
921 }
922 // Use throw directly to prevent file and line
923 throw TimeSteppingBreakdown{msg};
924 }
925}
926
927template<class TypeTag>
928template<class Solver>
929void
930AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
931checkTimeStepMinLimit_(const double new_time_step) const
932{
933 using Meas = UnitSystem::measure;
934 // If we have restarted (i.e. cut the timestep) too
935 // much, we have failed and throw an exception.
936 if (new_time_step < minTimeStep_()) {
937 std::string msg = "Solver failed to converge after cutting timestep to ";
938 if (Parameters::Get<Parameters::EnableTuning>()) {
939 const UnitSystem& unit_system = solver_().model().simulator().vanguard().eclState().getDeckUnitSystem();
940 msg += fmt::format(
941 "{:.3E} {}\nwhich is the minimum threshold given by the TUNING keyword\n",
942 unit_system.from_si(Meas::time, minTimeStep_()),
943 unit_system.name(Meas::time)
944 );
945 }
946 else {
947 msg += fmt::format(
948 "{:.3E} DAYS\nwhich is the minimum threshold given by option --solver-min-time-step\n",
949 minTimeStep_() / 86400.0
950 );
951 }
952 if (solverVerbose_()) {
953 OpmLog::error(msg);
954 }
955 // Use throw directly to prevent file and line
956 throw TimeSteppingBreakdown{msg};
957 }
958}
959
960template<class TypeTag>
961template<class Solver>
962void
963AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
964chopTimeStep_(const double new_time_step)
965{
966 setTimeStep_(new_time_step);
967 if (solverVerbose_()) {
968 const auto msg = fmt::format(fmt::runtime("{}\nTimestep chopped to {} days\n"),
969 this->cause_of_failure_,
970 unit::convert::to(this->substep_timer_.currentStepLength(), unit::day));
971 OpmLog::problem(msg);
972 }
973}
974
975template<class TypeTag>
976template<class Solver>
977bool
978AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
979chopTimeStepOrCloseFailingWells_(const double new_time_step)
980{
981 bool wells_shut = false;
982 // We are below the threshold, and will check if there are any
983 // wells that fails repeatedly (that means that it fails in the last three steps)
984 // we should close rather than chopping again.
985 // If we already have chopped the timestep two times that is
986 // new_time_step < minTimeStepBeforeClosingWells_()*restartFactor_()*restartFactor_()
987 // We also shut wells that fails only on this step.
988 const bool requireRepeatedFailures =
989 new_time_step > (minTimeStepBeforeClosingWells_() * restartFactor_() * restartFactor_());
990 const std::set<std::string> failing_wells =
991 detail::consistentlyFailingWells(solver_().model().stepReports(), requireRepeatedFailures);
992
993 if (failing_wells.empty()) {
994 // Found no wells to close, chop the timestep
995 chopTimeStep_(new_time_step);
996 } else {
997 // Close all consistently failing wells that are not under group control
998 std::vector<std::string> shut_wells;
999 for (const auto& well : failing_wells) {
1000 const bool was_shut =
1001 solver_().model().wellModel().forceShutWellByName(well,
1002 this->substep_timer_.simulationTimeElapsed(),
1003 /*dont_shut_grup_wells =*/ true);
1004 if (was_shut) {
1005 shut_wells.push_back(well);
1006 }
1007 }
1008 // If no wells are closed we also try to shut wells under group control
1009 if (shut_wells.empty()) {
1010 for (const auto& well : failing_wells) {
1011 const bool was_shut =
1012 solver_().model().wellModel().forceShutWellByName(well,
1013 this->substep_timer_.simulationTimeElapsed(),
1014 /*dont_shut_grup_wells =*/ false);
1015 if (was_shut) {
1016 shut_wells.push_back(well);
1017 }
1018 }
1019 }
1020 // If still no wells are closed we must fall back to chopping again
1021 if (shut_wells.empty()) {
1022 chopTimeStep_(new_time_step);
1023 } else {
1024 wells_shut = true;
1025 if (solverVerbose_()) {
1026 const std::string msg =
1027 fmt::format(fmt::runtime("\nProblematic well(s) were shut: {}"
1028 "(retrying timestep)\n"),
1029 fmt::join(shut_wells, " "));
1030 OpmLog::problem(msg);
1031 }
1032 }
1033 }
1034 return wells_shut;
1035}
1036
1037template<class TypeTag>
1038template<class Solver>
1039boost::posix_time::ptime
1040AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1041currentDateTime_() const
1042{
1043 return simulatorTimer_().currentDateTime();
1044}
1045
1046template<class TypeTag>
1047template<class Solver>
1048int
1049AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1050getNumIterations_(const SimulatorReportSingle &substep_report) const
1051{
1052 if (useNewtonIteration_()) {
1053 return substep_report.total_newton_iterations;
1054 }
1055 else {
1056 return substep_report.total_linear_iterations;
1057 }
1058}
1059
1060template<class TypeTag>
1061template<class Solver>
1062double
1063AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1064growthFactor_() const
1065{
1066 return this->adaptive_time_stepping_.growth_factor_;
1067}
1068
1069template<class TypeTag>
1070template<class Solver>
1071bool
1072AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1073ignoreConvergenceFailure_() const
1074{
1075 return adaptive_time_stepping_.ignore_convergence_failure_;
1076}
1077
1078template<class TypeTag>
1079template<class Solver>
1080bool
1081AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1082isReservoirCouplingMaster_() const
1083{
1084 return this->substepper_.isReservoirCouplingMaster_();
1085}
1086
1087template<class TypeTag>
1088template<class Solver>
1089bool
1090AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1091isReservoirCouplingSlave_() const
1092{
1093 return this->substepper_.isReservoirCouplingSlave_();
1094}
1095
1096template<class TypeTag>
1097template<class Solver>
1098void
1099AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1100markFirstSubStepAsFinished_() const
1101{
1102#ifdef RESERVOIR_COUPLING_ENABLED
1103 // Clear the first-substep flag after the first runSubStep_() call.
1104 // This ensures that master-slave synchronization only happens once per sync timestep,
1105 // not on retry attempts after convergence-driven timestep chops.
1106 if (isReservoirCouplingMaster_()) {
1107 reservoirCouplingMaster_().setFirstSubstepOfSyncTimestep(false);
1108 }
1109 else if (isReservoirCouplingSlave_()) {
1110 reservoirCouplingSlave_().setFirstSubstepOfSyncTimestep(false);
1111 }
1112#endif
1113 return;
1114}
1115
1116template<class TypeTag>
1117template<class Solver>
1118double
1119AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1120maxGrowth_() const
1121{
1122 return this->adaptive_time_stepping_.max_growth_;
1123}
1124
1125template<class TypeTag>
1126template<class Solver>
1127void
1128AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1129maybeReportSubStep_(SimulatorReportSingle substep_report) const
1130{
1131 if (timeStepVerbose_()) {
1132 std::ostringstream ss;
1133 substep_report.reportStep(ss);
1134 OpmLog::info(ss.str());
1135 }
1136}
1137
1138template<class TypeTag>
1139template<class Solver>
1140double
1141AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1142maybeRestrictTimeStepGrowth_(const double dt, double dt_estimate, const int restarts) const
1143{
1144 // limit the growth of the timestep size by the growth factor
1145 dt_estimate = std::min(dt_estimate, double(maxGrowth_() * dt));
1146 assert(dt_estimate > 0);
1147 // further restrict time step size growth after convergence problems
1148 if (restarts > 0) {
1149 dt_estimate = std::min(growthFactor_() * dt, dt_estimate);
1150 }
1151
1152 return dt_estimate;
1153}
1154
1155
1156template<class TypeTag>
1157template<class Solver>
1158void
1159AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1160maybeUpdateLastSubstepOfSyncTimestep_([[maybe_unused]] const double dt)
1161{
1162#ifdef RESERVOIR_COUPLING_ENABLED
1163 // For reservoir coupling slaves: predict if this substep will complete
1164 // the sync timestep. If so, timeStepSucceeded() will send production
1165 // data to the master (which is blocking on receive after its first substep).
1166 // This is used for summary data synchronization between slaves and master.
1167 if (isReservoirCouplingSlave_()) {
1168 const bool is_last = ReservoirCoupling::Seconds::compare_gt_or_eq(
1169 this->substep_timer_.simulationTimeElapsed() + dt,
1170 this->substep_timer_.totalTime()
1171 );
1172 reservoirCouplingSlave_().setLastSubstepOfSyncTimestep(is_last);
1173 }
1174#endif
1175}
1176
1177// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()
1178// It has to be called for each substep since TUNING might have been changed for next sub step due
1179// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
1180template<class TypeTag>
1181template<class Solver>
1182void
1183AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1184maybeUpdateTuningAndTimeStep_()
1185{
1186 // TODO: This function is currently only called if NEXTSTEP is activated from ACTIONX or
1187 // if the WCYCLE keyword needs to modify the current timestep. So this method should rather
1188 // be named maybeUpdateTimeStep_() or similar, since it should not update the tuning. However,
1189 // the current definition of the maybeUpdateTuning_() callback is actually calling
1190 // adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning) which is updating the tuning
1191 // see SimulatorFullyImplicitBlackoil::runStep() for more details.
1192 const auto old_value = suggestedNextTimestep_();
1193 if (this->substepper_.maybeUpdateTuning_(this->substep_timer_.simulationTimeElapsed(),
1194 this->substep_timer_.currentStepLength(),
1195 this->substep_timer_.currentStepNum()))
1196 {
1197 // Either NEXTSTEP and WCYCLE wants to change the current time step, but they cannot
1198 // change the current time step directly. Instead, they change the suggested next time step
1199 // by calling updateNEXTSTEP() via the maybeUpdateTuning() callback. We now need to update
1200 // the current time step to the new suggested time step and reset the suggested time step
1201 // to the old value.
1202 setTimeStep_(suggestedNextTimestep_());
1203 setSuggestedNextStep_(old_value);
1204 }
1205}
1206
1207template<class TypeTag>
1208template<class Solver>
1209double
1210AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1211minTimeStepBeforeClosingWells_() const
1212{
1213 return this->adaptive_time_stepping_.min_time_step_before_shutting_problematic_wells_;
1214}
1215
1216template<class TypeTag>
1217template<class Solver>
1218double
1219AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1220minTimeStep_() const
1221{
1222 return this->adaptive_time_stepping_.min_time_step_;
1223}
1224
1225#ifdef RESERVOIR_COUPLING_ENABLED
1226template<class TypeTag>
1227template<class Solver>
1228ReservoirCouplingMaster<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
1229AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1230reservoirCouplingMaster_() const
1231{
1232 return this->substepper_.reservoirCouplingMaster_();
1233}
1234
1235template<class TypeTag>
1236template<class Solver>
1237ReservoirCouplingSlave<typename AdaptiveTimeStepping<TypeTag>::Scalar>&
1238AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1239reservoirCouplingSlave_() const
1240{
1241 return this->substepper_.reservoirCouplingSlave_();
1242}
1243#endif
1244
1245template<class TypeTag>
1246template<class Solver>
1247double
1248AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1249restartFactor_() const
1250{
1251 return this->adaptive_time_stepping_.restart_factor_;
1252}
1253
1254template<class TypeTag>
1255template<class Solver>
1256SimulatorReportSingle
1257AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1258runSubStep_()
1259{
1260 SimulatorReportSingle substep_report;
1261
1262 auto handleFailure = [this, &substep_report]
1263 (const std::string& failure_reason, const std::exception& e, bool log_exception = true)
1264 {
1265 substep_report = solver_().failureReport();
1266 this->cause_of_failure_ = failure_reason;
1267 if (log_exception && solverVerbose_()) {
1268 OpmLog::debug(std::string("Caught Exception: ") + e.what());
1269 }
1270 };
1271
1272 try {
1273 substep_report = solver_().step(this->substep_timer_, &this->adaptive_time_stepping_.timeStepControl());
1274 if (solverVerbose_()) {
1275 // report number of linear iterations
1276 OpmLog::debug("Overall linear iterations used: "
1277 + std::to_string(substep_report.total_linear_iterations));
1278 }
1279 }
1280 catch (const TooManyIterations& e) {
1281 handleFailure("Solver convergence failure - Iteration limit reached", e);
1282 }
1283 catch (const TimeSteppingBreakdown& e) {
1284 handleFailure(e.what(), e);
1285 }
1286 catch (const ConvergenceMonitorFailure& e) {
1287 handleFailure("Convergence monitor failure", e, /*log_exception=*/false);
1288 }
1289 catch (const LinearSolverProblem& e) {
1290 handleFailure("Linear solver convergence failure", e);
1291 }
1292 catch (const NumericalProblem& e) {
1293 handleFailure("Solver convergence failure - Numerical problem encountered", e);
1294 }
1295 catch (const std::runtime_error& e) {
1296 handleFailure("Runtime error encountered", e);
1297 }
1298 catch (const Dune::ISTLError& e) {
1299 handleFailure("ISTL error - Time step too large", e);
1300 }
1301 catch (const Dune::MatrixBlockError& e) {
1302 handleFailure("Matrix block error", e);
1303 }
1304
1305 return substep_report;
1306}
1307
1308template<class TypeTag>
1309template<class Solver>
1310void
1311AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1312setTimeStep_(double dt_estimate)
1313{
1314 this->substep_timer_.provideTimeStepEstimate(dt_estimate);
1315}
1316
1317template<class TypeTag>
1318template<class Solver>
1319Solver&
1320AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1321solver_() const
1322{
1323 return this->substepper_.solver_;
1324}
1325
1326
1327template<class TypeTag>
1328template<class Solver>
1329int
1330AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1331solverRestartMax_() const
1332{
1333 return this->adaptive_time_stepping_.solver_restart_max_;
1334}
1335
1336template<class TypeTag>
1337template<class Solver>
1338void
1339AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1340setSuggestedNextStep_(double step)
1341{
1342 this->adaptive_time_stepping_.setSuggestedNextStep(step);
1343}
1344
1345template <class TypeTag>
1346template <class Solver>
1347const SimulatorTimer&
1348AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1349simulatorTimer_() const
1350{
1351 return this->substepper_.simulator_timer_;
1352}
1353
1354template <class TypeTag>
1355template <class Solver>
1356bool
1357AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1358solverVerbose_() const
1359{
1360 return this->adaptive_time_stepping_.solver_verbose_;
1361}
1362
1363template<class TypeTag>
1364template<class Solver>
1365boost::posix_time::ptime
1366AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1367startDateTime_() const
1368{
1369 return simulatorTimer_().startDateTime();
1370}
1371
1372template <class TypeTag>
1373template <class Solver>
1374double
1375AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1376suggestedNextTimestep_() const
1377{
1378 return this->adaptive_time_stepping_.suggestedNextStep();
1379}
1380
1381template <class TypeTag>
1382template <class Solver>
1383double
1384AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1385timeStepControlComputeEstimate_(const double dt, const int iterations,
1386 const AdaptiveSimulatorTimer& substepTimer) const
1387{
1388 // create object to compute the time error, simply forwards the call to the model
1389 const SolutionTimeErrorSolverWrapper<Solver> relative_change{solver_()};
1390 return this->adaptive_time_stepping_.time_step_control_->computeTimeStepSize(
1391 dt, iterations, relative_change, substepTimer);
1392}
1393
1394template <class TypeTag>
1395template <class Solver>
1396bool
1397AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1398timeStepVerbose_() const
1399{
1400 return this->adaptive_time_stepping_.timestep_verbose_;
1401}
1402
1403// The suggested time step is the stepsize that will be used as a first try for
1404// the next sub step. It is updated at the end of each substep. It can also be
1405// updated by the TUNING or NEXTSTEP keywords at the beginning of each report step or
1406// at the beginning of each substep by the ACTIONX keyword (via NEXTSTEP), this is
1407// done by the maybeUpdateTuning_() method which is called at the beginning of each substep
1408// (and the begginning of each report step). Note that the WCYCLE keyword can also update the
1409// suggested time step via the maybeUpdateTuning_() method.
1410template <class TypeTag>
1411template <class Solver>
1412void
1413AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1414updateSuggestedNextStep_()
1415{
1416 auto suggested_next_step = this->substep_timer_.currentStepLength();
1417 if (! std::isfinite(suggested_next_step)) { // check for NaN
1418 suggested_next_step = this->original_time_step_;
1419 }
1420 if (timeStepVerbose_()) {
1421 std::ostringstream ss;
1422 this->substep_timer_.report(ss);
1423 ss << "Suggested next step size = "
1424 << unit::convert::to(suggested_next_step, unit::day) << " (days)" << std::endl;
1425 OpmLog::debug(ss.str());
1426 }
1427 setSuggestedNextStep_(suggested_next_step);
1428}
1429
1430template <class TypeTag>
1431template <class Solver>
1432bool
1433AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1434useNewtonIteration_() const
1435{
1436 return this->adaptive_time_stepping_.use_newton_iteration_;
1437}
1438
1439template <class TypeTag>
1440template <class Solver>
1441double
1442AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1443writeOutput_() const
1444{
1445 time::StopWatch perf_timer;
1446 perf_timer.start();
1447 auto& problem = solver_().model().simulator().problem();
1448 problem.writeOutput(true);
1449 return perf_timer.secsSinceStart();
1450}
1451
1452/************************************************
1453 * Private class SolutionTimeErrorSolverWrapper
1454 * **********************************************/
1455
1456template<class TypeTag>
1457template<class Solver>
1458AdaptiveTimeStepping<TypeTag>::
1459SolutionTimeErrorSolverWrapper<Solver>::
1460SolutionTimeErrorSolverWrapper(const Solver& solver)
1461 : solver_{solver}
1462{}
1463
1464template<class TypeTag>
1465template<class Solver>
1466double AdaptiveTimeStepping<TypeTag>::SolutionTimeErrorSolverWrapper<Solver>::relativeChange() const
1467{
1468 // returns: || u^n+1 - u^n || / || u^n+1 ||
1469 return solver_.model().relativeChange();
1470}
1471
1472} // namespace Opm
1473
1474#endif // OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
Opm::AdaptiveTimeStepping
Adaptive time-stepping coordinator for the black-oil simulator.
Definition: AdaptiveTimeStepping.hpp:93
Opm::AdaptiveTimeStepping::max_growth_
double max_growth_
factor that limits the maximum growth of a time step
Definition: AdaptiveTimeStepping.hpp:422
Opm::AdaptiveTimeStepping::max_time_step_
double max_time_step_
maximal allowed time step size in days
Definition: AdaptiveTimeStepping.hpp:423
Opm::AdaptiveTimeStepping::solver_verbose_
bool solver_verbose_
solver verbosity
Definition: AdaptiveTimeStepping.hpp:427
Opm::AdaptiveTimeStepping::solver_restart_max_
int solver_restart_max_
how many restart of solver are allowed
Definition: AdaptiveTimeStepping.hpp:426
Opm::AdaptiveTimeStepping::timestep_after_event_
double timestep_after_event_
suggested size of timestep after an event
Definition: AdaptiveTimeStepping.hpp:431
Opm::AdaptiveTimeStepping::init_
void init_(const UnitSystem &unitSystem)
Definition: AdaptiveTimeStepping_impl.hpp:428
Opm::AdaptiveTimeStepping::setSuggestedNextStep
void setSuggestedNextStep(const double x)
Set the suggested length for the next substep [s].
Definition: AdaptiveTimeStepping_impl.hpp:299
Opm::AdaptiveTimeStepping::suggestedNextStep
double suggestedNextStep() const
Definition: AdaptiveTimeStepping_impl.hpp:307
Opm::AdaptiveTimeStepping::operator==
bool operator==(const AdaptiveTimeStepping< TypeTag > &rhs) const
Definition: AdaptiveTimeStepping_impl.hpp:139
Opm::AdaptiveTimeStepping::serializationTestObjectSimple
static AdaptiveTimeStepping< TypeTag > serializationTestObjectSimple()
Definition: AdaptiveTimeStepping_impl.hpp:282
Opm::AdaptiveTimeStepping::ignore_convergence_failure_
bool ignore_convergence_failure_
continue instead of stop when minimum time step is reached
Definition: AdaptiveTimeStepping.hpp:425
Opm::AdaptiveTimeStepping::serializeOp
void serializeOp(Serializer &serializer)
Definition: AdaptiveTimeStepping_impl.hpp:211
Opm::AdaptiveTimeStepping::updateTUNING
void updateTUNING(double max_next_tstep, const Tuning &tuning)
Apply TUNING keyword parameters.
Definition: AdaptiveTimeStepping_impl.hpp:336
Opm::AdaptiveTimeStepping::time_step_control_type_
TimeStepControlType time_step_control_type_
type of time step control object
Definition: AdaptiveTimeStepping.hpp:418
Opm::AdaptiveTimeStepping::timeStepControl
const TimeStepControlInterface & timeStepControl() const
Definition: AdaptiveTimeStepping_impl.hpp:315
Opm::AdaptiveTimeStepping::TuningUpdateCallback
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
Opm::AdaptiveTimeStepping::full_timestep_initially_
bool full_timestep_initially_
beginning with the size of the time step from data file
Definition: AdaptiveTimeStepping.hpp:430
Opm::AdaptiveTimeStepping::step
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:196
Opm::AdaptiveTimeStepping::growth_factor_
double growth_factor_
factor to multiply time step when solver recovered from failed convergence
Definition: AdaptiveTimeStepping.hpp:421
Opm::AdaptiveTimeStepping::restart_factor_
double restart_factor_
factor to multiply time step with when solver fails to converge
Definition: AdaptiveTimeStepping.hpp:420
Opm::AdaptiveTimeStepping::min_time_step_
double min_time_step_
minimal allowed time step size before throwing
Definition: AdaptiveTimeStepping.hpp:424
Opm::AdaptiveTimeStepping::updateNEXTSTEP
void updateNEXTSTEP(double max_next_tstep)
Set suggested_next_timestep_ to max_next_tstep iff max_next_tstep > 0.
Definition: AdaptiveTimeStepping_impl.hpp:324
Opm::AdaptiveTimeStepping::serializationTestObjectHardcoded
static AdaptiveTimeStepping< TypeTag > serializationTestObjectHardcoded()
Definition: AdaptiveTimeStepping_impl.hpp:258
Opm::AdaptiveTimeStepping::time_step_control_
TimeStepController time_step_control_
time step control object
Definition: AdaptiveTimeStepping.hpp:419
Opm::AdaptiveTimeStepping::serializationTestObjectPIDIt
static AdaptiveTimeStepping< TypeTag > serializationTestObjectPIDIt()
Definition: AdaptiveTimeStepping_impl.hpp:274
Opm::AdaptiveTimeStepping::min_time_step_before_shutting_problematic_wells_
double min_time_step_before_shutting_problematic_wells_
< shut problematic wells when time step size in days are less than this
Definition: AdaptiveTimeStepping.hpp:435
Opm::AdaptiveTimeStepping::serializationTestObject3rdOrder
static AdaptiveTimeStepping< TypeTag > serializationTestObject3rdOrder()
Definition: AdaptiveTimeStepping_impl.hpp:290
Opm::AdaptiveTimeStepping::report
SimulatorReport & report()
Definition: AdaptiveTimeStepping_impl.hpp:250
Opm::AdaptiveTimeStepping::serializationTestObjectPID
static AdaptiveTimeStepping< TypeTag > serializationTestObjectPID()
Definition: AdaptiveTimeStepping_impl.hpp:266
Opm::AdaptiveTimeStepping::registerParameters
static void registerParameters()
Definition: AdaptiveTimeStepping_impl.hpp:185
Opm::AdaptiveTimeStepping::use_newton_iteration_
bool use_newton_iteration_
use newton iteration count for adaptive time step control
Definition: AdaptiveTimeStepping.hpp:432
Opm::SimulatorTimer
Definition: SimulatorTimer.hpp:39
Opm::TimeStepControlInterface
Definition: TimeStepControlInterface.hpp:51
Opm::Parameters::Get
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition: parametersystem.hpp:187
Opm::detail::logTimer
void logTimer(const AdaptiveSimulatorTimer &substep_timer)
Opm::detail::registerAdaptiveParameters
void registerAdaptiveParameters()
Opm::detail::consistentlyFailingWells
std::set< std::string > consistentlyFailingWells(const std::vector< StepReport > &sr, bool requireRepeatedFailures)
Opm::detail::createController
std::tuple< TimeStepControlType, std::unique_ptr< TimeStepControlInterface >, bool > createController(const UnitSystem &unitSystem)
Opm
Definition: blackoilbioeffectsmodules.hh:45
Opm::to_string
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
parametersystem.hpp
This file provides the infrastructure to retrieve run-time parameters.
Opm::ReservoirCoupling::Seconds::compare_gt_or_eq
static bool compare_gt_or_eq(double a, double b)
Determines if a is greater than b within the specified tolerance.
Opm::SimulatorReport
Definition: SimulatorReport.hpp:122