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