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#ifdef RESERVOIR_COUPLING_ENABLED
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#endif
529
530template<class TypeTag>
531template<class Solver>
532void
533AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
534maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step)
535{
536 this->adaptive_time_stepping_.maybeModifySuggestedTimeStepAtBeginningOfReportStep_(
537 original_time_step, this->is_event_
538 );
539}
540
541// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()
542// It has to be called for each substep since TUNING might have been changed for next sub step due
543// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
544template<class TypeTag>
545template<class Solver>
546bool
547AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
548maybeUpdateTuning_(double elapsed, double dt, int sub_step_number) const
549{
550 return this->tuning_updater_(elapsed, dt, sub_step_number);
551}
552
553template<class TypeTag>
554template<class Solver>
555double
556AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
557maxTimeStep_() const
558{
559 return this->adaptive_time_stepping_.max_time_step_;
560}
561
562template <class TypeTag>
563template <class Solver>
564SimulatorReport
565AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
566runStepOriginal_()
567{
568 auto elapsed = this->simulator_timer_.simulationTimeElapsed();
569 auto original_time_step = this->simulator_timer_.currentStepLength();
570 auto report_step = this->simulator_timer_.reportStepNum();
571 maybeUpdateTuning_(elapsed, original_time_step, report_step);
572 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(original_time_step);
573
574 AdaptiveSimulatorTimer substep_timer{
575 this->simulator_timer_.startDateTime(),
577 elapsed,
578 suggestedNextTimestep_(),
579 report_step,
580 maxTimeStep_()
581 };
582 SubStepIteration<Solver> substepIteration{*this, substep_timer, original_time_step, /*final_step=*/true};
583 return substepIteration.run();
584}
585
586#ifdef RESERVOIR_COUPLING_ENABLED
587template <class TypeTag>
588template <class Solver>
589ReservoirCouplingMaster&
590AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
591reservoirCouplingMaster_()
592{
593 return *adaptive_time_stepping_.reservoir_coupling_master_;
594}
595#endif
596
597#ifdef RESERVOIR_COUPLING_ENABLED
598template <class TypeTag>
599template <class Solver>
600ReservoirCouplingSlave&
601AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
602reservoirCouplingSlave_()
603{
604 return *this->adaptive_time_stepping_.reservoir_coupling_slave_;
605}
606#endif
607
608#ifdef RESERVOIR_COUPLING_ENABLED
609// Description of the reservoir coupling master and slave substep loop
610// -------------------------------------------------------------------
611// The master and slave processes attempts to reach the end of the report step using a series of substeps
612// (also called timesteps). Each substep have an upper limit that is roughly determined by a combination
613// of the keywords TUNING (through the TSINIT, TSMAXZ values), NEXSTEP, WCYCLE, and the start of the
614// next report step (the last substep needs to coincide with this time). Note that NEXTSTEP can be
615// updated from an ACTIONX keyword. Although, this comment focuses on the maximum substep limit,
616// note that there is also a lower limit on the substep length. And the substep sizes will be adjusted
617// automatically (or retried) based on the convergence behavior of the solver and other criteria.
618//
619// When using reservoir coupling, the upper limit on each substep is further limited to: a) not overshoot
620// next report date of a slave reservoir, and b) to keep the flow rate of the slave groups within
621// certain limits. To determine this potential further limitation on the substep, the following procedure
622// is used at the beginning of each master substep:
623// - First, the slaves sends the master the date of their next report step
624// - The master receives the date of the next report step from the slaves
625// - If necessary, the master computes a modified substep length based on the received dates
626// TODO: explain how the substep is limited to keep the flow rate of the slave groups within certain
627// limits. Since this is not implemented yet, this part of the procedure is not explained here.
628//
629// Then, after the master has determined the substep length using the above procedure, it sends it
630// to the slaves. The slaves will now use the end of this substep as a fixed point (like a mini report
631// step), that they will try to reach using a single substep or multiple substeps. The end of the
632// substep received from the master will either conincide with the end of its current report step, or
633// it will end before (it cannot not end after since the master has already adjusted the step to not
634// exceed any report time in a slave). If the end of this substep is earlier in time than its next report
635// date, the slave will enter a loop and wait for the master to send it a new substep.
636// The loop is finished when the end of the received time step conincides with the end of its current
637// report step.
638
639template <class TypeTag>
640template <class Solver>
641SimulatorReport
642AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
643runStepReservoirCouplingMaster_()
644{
645 bool substep_done = false;
646 int iteration = 0;
647 const double original_time_step = this->simulator_timer_.currentStepLength();
648 double current_time{this->simulator_timer_.simulationTimeElapsed()};
651 SimulatorReport report;
652 while(!substep_done) {
653 reservoirCouplingMaster_().receiveNextReportDateFromSlaves();
654 if (iteration == 0) {
655 maybeUpdateTuning_(current_time, current_step_length, /*substep=*/0);
656 }
659 reservoirCouplingMaster_().sendNextTimeStepToSlaves(current_step_length);
660 if (iteration == 0) {
661 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(current_step_length);
662 }
663 AdaptiveSimulatorTimer substep_timer{
664 this->simulator_timer_.startDateTime(),
665 /*stepLength=*/current_step_length,
666 /*elapsedTime=*/current_time,
667 /*timeStepEstimate=*/suggestedNextTimestep_(),
668 /*reportStep=*/this->simulator_timer_.reportStepNum(),
669 maxTimeStep_()
670 };
673 );
674 SubStepIteration<Solver> substepIteration{*this, substep_timer, current_step_length, final_step};
676 report += sub_steps_report;
678 if (final_step) {
679 break;
680 }
681 iteration++;
682 }
683 return report;
684}
685#endif
686
687#ifdef RESERVOIR_COUPLING_ENABLED
688template <class TypeTag>
689template <class Solver>
690SimulatorReport
691AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
692runStepReservoirCouplingSlave_()
693{
694 bool substep_done = false;
695 int iteration = 0;
696 const double original_time_step = this->simulator_timer_.currentStepLength();
697 double current_time{this->simulator_timer_.simulationTimeElapsed()};
699 SimulatorReport report;
700 while(!substep_done) {
701 reservoirCouplingSlave_().sendNextReportDateToMasterProcess();
702 auto timestep = reservoirCouplingSlave_().receiveNextTimeStepFromMaster();
703 if (iteration == 0) {
704 maybeUpdateTuning_(current_time, original_time_step, /*substep=*/0);
705 maybeModifySuggestedTimeStepAtBeginningOfReportStep_(timestep);
706 }
707 AdaptiveSimulatorTimer substep_timer{
708 this->simulator_timer_.startDateTime(),
709 /*step_length=*/timestep,
710 /*elapsed_time=*/current_time,
711 /*time_step_estimate=*/suggestedNextTimestep_(),
712 this->simulator_timer_.reportStepNum(),
713 maxTimeStep_()
714 };
717 );
718 SubStepIteration<Solver> substepIteration{*this, substep_timer, timestep, final_step};
720 report += sub_steps_report;
722 if (final_step) {
723 substep_done = true;
724 break;
725 }
726 iteration++;
727 }
728 return report;
729}
730
731#endif
732
733template <class TypeTag>
734template <class Solver>
735double
736AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
737suggestedNextTimestep_() const
738{
739 return this->adaptive_time_stepping_.suggestedNextStep();
740}
741
742
743
744/************************************************
745 * Private class SubStepIteration public methods
746 ************************************************/
747
748template<class TypeTag>
749template<class Solver>
750AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
751SubStepIteration(
753 AdaptiveSimulatorTimer& substep_timer,
754 const double original_time_step,
755 bool final_step
756)
757 : substepper_{substepper}
758 , substep_timer_{substep_timer}
759 , original_time_step_{original_time_step}
760 , final_step_{final_step}
761 , adaptive_time_stepping_{substepper.getAdaptiveTimerStepper()}
762{
763}
764
765template <class TypeTag>
766template <class Solver>
767SimulatorReport
768AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
769run()
770{
771 auto& simulator = solver_().model().simulator();
772 auto& problem = simulator.problem();
773 // counter for solver restarts
774 int restarts = 0;
775 SimulatorReport report;
776
777 // sub step time loop
778 while (!this->substep_timer_.done()) {
779 // if we just chopped the timestep due to convergence i.e. restarts>0
780 // we dont what to update the next timestep based on Tuning
781 if (restarts == 0) {
782 maybeUpdateTuningAndTimeStep_();
783 }
784 const double dt = this->substep_timer_.currentStepLength();
785 if (timeStepVerbose_()) {
786 detail::logTimer(this->substep_timer_);
787 }
788
789 auto substep_report = runSubStep_();
790
791 //Pass substep to eclwriter for summary output
792 problem.setSubStepReport(substep_report);
793 auto& full_report = adaptive_time_stepping_.report();
795 problem.setSimulationReport(full_report);
796
797 report += substep_report;
798
799 if (substep_report.converged || checkContinueOnUnconvergedSolution_(dt)) {
800 ++this->substep_timer_; // advance by current dt
801
802 const int iterations = getNumIterations_(substep_report);
803 auto dt_estimate = timeStepControlComputeEstimate_(
804 dt, iterations, this->substep_timer_);
805
806 assert(dt_estimate > 0);
807 dt_estimate = maybeRestrictTimeStepGrowth_(dt, dt_estimate, restarts);
808 restarts = 0; // solver converged, reset restarts counter
809
810 maybeReportSubStep_(substep_report);
811 if (this->final_step_ && this->substep_timer_.done()) {
812 // if the time step is done we do not need to write it as this will be done
813 // by the simulator anyway.
814 }
815 else {
816 report.success.output_write_time += writeOutput_();
817 }
818
819 // set new time step length
820 setTimeStep_(dt_estimate);
821
822 report.success.converged = this->substep_timer_.done();
823 this->substep_timer_.setLastStepFailed(false);
824 }
825 else { // in case of no convergence
826 this->substep_timer_.setLastStepFailed(true);
827 checkTimeStepMaxRestartLimit_(restarts);
828
829 const double new_time_step = restartFactor_() * dt;
830 checkTimeStepMinLimit_(new_time_step);
831 bool wells_shut = false;
832 if (new_time_step > minTimeStepBeforeClosingWells_()) {
833 chopTimeStep_(new_time_step);
834 } else {
835 wells_shut = chopTimeStepOrCloseFailingWells_(new_time_step);
836 }
837 if (wells_shut) {
838 setTimeStep_(dt); // retry the old timestep
839 }
840 else {
841 restarts++; // only increase if no wells were shut
842 }
843 }
844 problem.setNextTimeStepSize(this->substep_timer_.currentStepLength());
845 }
846 updateSuggestedNextStep_();
847 return report;
848}
849
850
851/************************************************
852 * Private class SubStepIteration private methods
853 ************************************************/
854
855
856template<class TypeTag>
857template<class Solver>
858bool
859AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
860checkContinueOnUnconvergedSolution_(double dt) const
861{
862 bool continue_on_uncoverged_solution = ignoreConvergenceFailure_() && dt <= minTimeStep_();
863 if (continue_on_uncoverged_solution && solverVerbose_()) {
864 // NOTE: This method is only called if the solver failed to converge.
865 const auto msg = fmt::format(
866 "Solver failed to converge but timestep {} is smaller or equal to {}\n"
867 "which is the minimum threshold given by option --solver-min-time-step\n",
868 dt, minTimeStep_()
869 );
870 OpmLog::problem(msg);
871 }
873}
874
875template<class TypeTag>
876template<class Solver>
877void
878AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
879checkTimeStepMaxRestartLimit_(const int restarts) const
880{
881 // If we have restarted (i.e. cut the timestep) too
882 // many times, we have failed and throw an exception.
883 if (restarts >= solverRestartMax_()) {
884 const auto msg = fmt::format(
885 "Solver failed to converge after cutting timestep {} times.", restarts
886 );
887 if (solverVerbose_()) {
888 OpmLog::error(msg);
889 }
890 // Use throw directly to prevent file and line
892 }
893}
894
895template<class TypeTag>
896template<class Solver>
897void
898AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
899checkTimeStepMinLimit_(const double new_time_step) const
900{
901 // If we have restarted (i.e. cut the timestep) too
902 // much, we have failed and throw an exception.
903 if (new_time_step < minTimeStep_()) {
904 const auto msg = fmt::format(
905 "Solver failed to converge after cutting timestep to {}\n"
906 "which is the minimum threshold given by option --solver-min-time-step\n",
907 minTimeStep_()
908 );
909 if (solverVerbose_()) {
910 OpmLog::error(msg);
911 }
912 // Use throw directly to prevent file and line
914 }
915}
916
917template<class TypeTag>
918template<class Solver>
919void
920AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
921chopTimeStep_(const double new_time_step)
922{
923 setTimeStep_(new_time_step);
924 if (solverVerbose_()) {
925 const auto msg = fmt::format("{}\nTimestep chopped to {} days\n",
926 this->cause_of_failure_,
927 unit::convert::to(this->substep_timer_.currentStepLength(), unit::day));
928 OpmLog::problem(msg);
929 }
930}
931
932template<class TypeTag>
933template<class Solver>
934bool
935AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
936chopTimeStepOrCloseFailingWells_(const double new_time_step)
937{
938 bool wells_shut = false;
939 // We are below the threshold, and will check if there are any
940 // wells that fails repeatedly (that means that it fails in the last three steps)
941 // we should close rather than chopping again.
942 // If we already have chopped the timestep two times that is
943 // new_time_step < minTimeStepBeforeClosingWells_()*restartFactor_()*restartFactor_()
944 // We also shut wells that fails only on this step.
945 bool requireRepeatedFailures = new_time_step > ( minTimeStepBeforeClosingWells_()*restartFactor_()*restartFactor_());
946 std::set<std::string> failing_wells = detail::consistentlyFailingWells(solver_().model().stepReports(), requireRepeatedFailures);
947
948 if (failing_wells.empty()) {
949 // Found no wells to close, chop the timestep
950 chopTimeStep_(new_time_step);
951 } else {
952 // Close all consistently failing wells that are not under group control
953 std::vector<std::string> shut_wells;
954 for (const auto& well : failing_wells) {
955 bool was_shut = solver_().model().wellModel().forceShutWellByName(
956 well, this->substep_timer_.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ true);
957 if (was_shut) {
958 shut_wells.push_back(well);
959 }
960 }
961 // If no wells are closed we also try to shut wells under group control
962 if (shut_wells.empty()) {
963 for (const auto& well : failing_wells) {
964 bool was_shut = solver_().model().wellModel().forceShutWellByName(
965 well, this->substep_timer_.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ false);
966 if (was_shut) {
967 shut_wells.push_back(well);
968 }
969 }
970 }
971 // If still no wells are closed we must fall back to chopping again
972 if (shut_wells.empty()) {
973 chopTimeStep_(new_time_step);
974 } else {
975 wells_shut = true;
976 if (solverVerbose_()) {
977 const std::string msg =
978 fmt::format("\nProblematic well(s) were shut: {}"
979 "(retrying timestep)\n",
980 fmt::join(shut_wells, " "));
981 OpmLog::problem(msg);
982 }
983 }
984 }
985 return wells_shut;
986}
987
988template<class TypeTag>
989template<class Solver>
990boost::posix_time::ptime
991AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
992currentDateTime_() const
993{
994 return simulatorTimer_().currentDateTime();
995}
996
997template<class TypeTag>
998template<class Solver>
999int
1000AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1001getNumIterations_(const SimulatorReportSingle &substep_report) const
1002{
1003 if (useNewtonIteration_()) {
1004 return substep_report.total_newton_iterations;
1005 }
1006 else {
1007 return substep_report.total_linear_iterations;
1008 }
1009}
1010
1011template<class TypeTag>
1012template<class Solver>
1013double
1014AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1015growthFactor_() const
1016{
1017 return this->adaptive_time_stepping_.growth_factor_;
1018}
1019
1020template<class TypeTag>
1021template<class Solver>
1022bool
1023AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1024ignoreConvergenceFailure_() const
1025{
1026 return adaptive_time_stepping_.ignore_convergence_failure_;
1027}
1028
1029template<class TypeTag>
1030template<class Solver>
1031double
1032AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1033maxGrowth_() const
1034{
1035 return this->adaptive_time_stepping_.max_growth_;
1036}
1037
1038template<class TypeTag>
1039template<class Solver>
1040void
1041AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1042maybeReportSubStep_(SimulatorReportSingle substep_report) const
1043{
1044 if (timeStepVerbose_()) {
1045 std::ostringstream ss;
1046 substep_report.reportStep(ss);
1047 OpmLog::info(ss.str());
1048 }
1049}
1050
1051template<class TypeTag>
1052template<class Solver>
1053double
1054AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1055maybeRestrictTimeStepGrowth_(const double dt, double dt_estimate, const int restarts) const
1056{
1057 // limit the growth of the timestep size by the growth factor
1058 dt_estimate = std::min(dt_estimate, double(maxGrowth_() * dt));
1059 assert(dt_estimate > 0);
1060 // further restrict time step size growth after convergence problems
1061 if (restarts > 0) {
1062 dt_estimate = std::min(growthFactor_() * dt, dt_estimate);
1063 }
1064
1065 return dt_estimate;
1066}
1067
1068// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()
1069// It has to be called for each substep since TUNING might have been changed for next sub step due
1070// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
1071template<class TypeTag>
1072template<class Solver>
1073void
1074AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1075maybeUpdateTuningAndTimeStep_()
1076{
1077 // TODO: This function is currently only called if NEXTSTEP is activated from ACTIONX or
1078 // if the WCYCLE keyword needs to modify the current timestep. So this method should rather
1079 // be named maybeUpdateTimeStep_() or similar, since it should not update the tuning. However,
1080 // the current definition of the maybeUpdateTuning_() callback is actually calling
1081 // adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning) which is updating the tuning
1082 // see SimulatorFullyImplicitBlackoil::runStep() for more details.
1083 auto old_value = suggestedNextTimestep_();
1084 if (this->substepper_.maybeUpdateTuning_(this->substep_timer_.simulationTimeElapsed(),
1085 this->substep_timer_.currentStepLength(),
1086 this->substep_timer_.currentStepNum()))
1087 {
1088 // Either NEXTSTEP and WCYCLE wants to change the current time step, but they cannot
1089 // change the current time step directly. Instead, they change the suggested next time step
1090 // by calling updateNEXTSTEP() via the maybeUpdateTuning() callback. We now need to update
1091 // the current time step to the new suggested time step and reset the suggested time step
1092 // to the old value.
1093 setTimeStep_(suggestedNextTimestep_());
1094 setSuggestedNextStep_(old_value);
1095 }
1096}
1097
1098template<class TypeTag>
1099template<class Solver>
1100double
1101AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1102minTimeStepBeforeClosingWells_() const
1103{
1104 return this->adaptive_time_stepping_.min_time_step_before_shutting_problematic_wells_;
1105}
1106
1107template<class TypeTag>
1108template<class Solver>
1109double
1110AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1111minTimeStep_() const
1112{
1113 return this->adaptive_time_stepping_.min_time_step_;
1114}
1115
1116template<class TypeTag>
1117template<class Solver>
1118double
1119AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1120restartFactor_() const
1121{
1122 return this->adaptive_time_stepping_.restart_factor_;
1123}
1124
1125template<class TypeTag>
1126template<class Solver>
1127SimulatorReportSingle
1128AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1129runSubStep_()
1130{
1131 SimulatorReportSingle substep_report;
1132
1133 auto handleFailure = [this, &substep_report]
1134 (const std::string& failure_reason, const std::exception& e, bool log_exception = true)
1135 {
1136 substep_report = solver_().failureReport();
1137 this->cause_of_failure_ = failure_reason;
1138 if (log_exception && solverVerbose_()) {
1139 std::string message;
1140 message = "Caught Exception: ";
1141 message += e.what();
1142 OpmLog::debug(message);
1143 }
1144 };
1145
1146 try {
1147 substep_report = solver_().step(this->substep_timer_);
1148 if (solverVerbose_()) {
1149 // report number of linear iterations
1150 OpmLog::debug("Overall linear iterations used: "
1151 + std::to_string(substep_report.total_linear_iterations));
1152 }
1153 }
1154 catch (const TooManyIterations& e) {
1155 handleFailure("Solver convergence failure - Iteration limit reached", e);
1156 }
1157 catch (const ConvergenceMonitorFailure& e) {
1158 handleFailure("Convergence monitor failure", e, /*log_exception=*/false);
1159 }
1160 catch (const LinearSolverProblem& e) {
1161 handleFailure("Linear solver convergence failure", e);
1162 }
1163 catch (const NumericalProblem& e) {
1164 handleFailure("Solver convergence failure - Numerical problem encountered", e);
1165 }
1166 catch (const std::runtime_error& e) {
1167 handleFailure("Runtime error encountered", e);
1168 }
1169 catch (const Dune::ISTLError& e) {
1170 handleFailure("ISTL error - Time step too large", e);
1171 }
1172 catch (const Dune::MatrixBlockError& e) {
1173 handleFailure("Matrix block error", e);
1174 }
1175
1176 return substep_report;
1177}
1178
1179template<class TypeTag>
1180template<class Solver>
1181void
1182AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1183setTimeStep_(double dt_estimate)
1184{
1185 this->substep_timer_.provideTimeStepEstimate(dt_estimate);
1186}
1187
1188template<class TypeTag>
1189template<class Solver>
1190Solver&
1191AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1192solver_() const
1193{
1194 return this->substepper_.solver_;
1195}
1196
1197
1198template<class TypeTag>
1199template<class Solver>
1200int
1201AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1202solverRestartMax_() const
1203{
1204 return this->adaptive_time_stepping_.solver_restart_max_;
1205}
1206
1207template<class TypeTag>
1208template<class Solver>
1209void
1210AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1211setSuggestedNextStep_(double step)
1212{
1213 this->adaptive_time_stepping_.setSuggestedNextStep(step);
1214}
1215
1216template <class TypeTag>
1217template <class Solver>
1218const SimulatorTimer&
1219AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1220simulatorTimer_() const
1221{
1222 return this->substepper_.simulator_timer_;
1223}
1224
1225template <class TypeTag>
1226template <class Solver>
1227bool
1228AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1229solverVerbose_() const
1230{
1231 return this->adaptive_time_stepping_.solver_verbose_;
1232}
1233
1234template<class TypeTag>
1235template<class Solver>
1236boost::posix_time::ptime
1237AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1238startDateTime_() const
1239{
1240 return simulatorTimer_().startDateTime();
1241}
1242
1243template <class TypeTag>
1244template <class Solver>
1245double
1246AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1247suggestedNextTimestep_() const
1248{
1249 return this->adaptive_time_stepping_.suggestedNextStep();
1250}
1251
1252template <class TypeTag>
1253template <class Solver>
1254double
1255AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1256timeStepControlComputeEstimate_(const double dt, const int iterations, AdaptiveSimulatorTimer& substepTimer) const
1257{
1258 // create object to compute the time error, simply forwards the call to the model
1260 return this->adaptive_time_stepping_.time_step_control_->computeTimeStepSize(
1261 dt, iterations, relative_change, substepTimer);
1262}
1263
1264template <class TypeTag>
1265template <class Solver>
1266bool
1267AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1268timeStepVerbose_() const
1269{
1270 return this->adaptive_time_stepping_.timestep_verbose_;
1271}
1272
1273// The suggested time step is the stepsize that will be used as a first try for
1274// the next sub step. It is updated at the end of each substep. It can also be
1275// updated by the TUNING or NEXTSTEP keywords at the beginning of each report step or
1276// at the beginning of each substep by the ACTIONX keyword (via NEXTSTEP), this is
1277// done by the maybeUpdateTuning_() method which is called at the beginning of each substep
1278// (and the begginning of each report step). Note that the WCYCLE keyword can also update the
1279// suggested time step via the maybeUpdateTuning_() method.
1280template <class TypeTag>
1281template <class Solver>
1282void
1283AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1284updateSuggestedNextStep_()
1285{
1286 auto suggested_next_step = this->substep_timer_.currentStepLength();
1287 if (! std::isfinite(suggested_next_step)) { // check for NaN
1288 suggested_next_step = this->original_time_step_;
1289 }
1290 if (timeStepVerbose_()) {
1291 std::ostringstream ss;
1292 this->substep_timer_.report(ss);
1293 ss << "Suggested next step size = "
1294 << unit::convert::to(suggested_next_step, unit::day) << " (days)" << std::endl;
1295 OpmLog::debug(ss.str());
1296 }
1297 setSuggestedNextStep_(suggested_next_step);
1298}
1299
1300template <class TypeTag>
1301template <class Solver>
1302bool
1303AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1304useNewtonIteration_() const
1305{
1306 return this->adaptive_time_stepping_.use_newton_iteration_;
1307}
1308
1309template <class TypeTag>
1310template <class Solver>
1311double
1312AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
1313writeOutput_() const
1314{
1315 time::StopWatch perf_timer;
1316 perf_timer.start();
1317 auto& problem = solver_().model().simulator().problem();
1318 problem.writeOutput(true);
1319 return perf_timer.secsSinceStart();
1320}
1321
1322/************************************************
1323 * Private class SolutionTimeErrorSolverWrapper
1324 * **********************************************/
1325
1326template<class TypeTag>
1327template<class Solver>
1328AdaptiveTimeStepping<TypeTag>::SolutionTimeErrorSolverWrapper<Solver>::SolutionTimeErrorSolverWrapper(
1329 const Solver& solver
1330)
1331 : solver_{solver}
1332{}
1333
1334template<class TypeTag>
1335template<class Solver>
1336double AdaptiveTimeStepping<TypeTag>::SolutionTimeErrorSolverWrapper<Solver>::relativeChange() const
1337{
1338 // returns: || u^n+1 - u^n || / || u^n+1 ||
1339 return solver_.model().relativeChange();
1340}
1341
1342} // namespace Opm
1343
1344#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:119