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