19 #ifndef OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
20 #define OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
26 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
27 #include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
28 #include <opm/simulators/timestepping/TimeStepControl.hpp>
29 #include <opm/core/utility/StopWatch.hpp>
30 #include <opm/common/Exceptions.hpp>
31 #include <opm/common/OpmLog/OpmLog.hpp>
32 #include <dune/istl/istlexception.hh>
33 #include <dune/istl/ilu.hh>
34 #include <opm/parser/eclipse/EclipseState/Schedule/Tuning.hpp>
40 template <
class Solver,
class State>
43 const Solver& solver_;
44 const State& previous_;
45 const State& current_;
48 const State& previous,
49 const State& current )
51 previous_( previous ),
58 return solver_.model().relativeChange( previous_, current_ );
63 void logException(
const E& exception,
bool verbose)
68 message =
"Caught Exception: ";
69 message += exception.what();
70 OpmLog::debug(message);
79 const ParameterGroup& param,
80 const bool terminal_output )
82 , restart_factor_( tuning.getTSFCNV(time_step) )
83 , growth_factor_(tuning.getTFDIFF(time_step) )
84 , max_growth_( tuning.getTSFMAX(time_step) )
86 , max_time_step_( tuning.getTSMAXZ(time_step) )
87 , solver_restart_max_( param.getDefault(
"solver.restart", int(10) ) )
88 , solver_verbose_( param.getDefault(
"solver.verbose", bool(true) ) && terminal_output )
89 , timestep_verbose_( param.getDefault(
"timestep.verbose", bool(true) ) && terminal_output )
90 , suggested_next_timestep_( tuning.getTSINIT(time_step) )
91 , full_timestep_initially_( param.getDefault(
"full_timestep_initially", bool(false) ) )
92 , timestep_after_event_( tuning.getTMAXWC(time_step))
93 , use_newton_iteration_(false)
100 const bool terminal_output )
102 , restart_factor_( param.getDefault(
"solver.restartfactor", double(0.33) ) )
103 , growth_factor_( param.getDefault(
"solver.growthfactor", double(2) ) )
104 , max_growth_( param.getDefault(
"timestep.control.maxgrowth", double(3.0) ) )
106 , max_time_step_( unit::convert::from(param.getDefault(
"timestep.max_timestep_in_days", 365.0 ), unit::day) )
107 , solver_restart_max_( param.getDefault(
"solver.restart", int(10) ) )
108 , solver_verbose_( param.getDefault(
"solver.verbose", bool(true) ) && terminal_output )
109 , timestep_verbose_( param.getDefault(
"timestep.verbose", bool(true) ) && terminal_output )
110 , suggested_next_timestep_( unit::convert::from(param.getDefault(
"timestep.initial_timestep_in_days", -1.0 ), unit::day) )
111 , full_timestep_initially_( param.getDefault(
"full_timestep_initially", bool(false) ) )
112 , timestep_after_event_( unit::convert::from(param.getDefault(
"timestep.timestep_in_days_after_event", -1.0 ), unit::day))
113 , use_newton_iteration_(false)
118 inline void AdaptiveTimeStepping::
119 init(
const ParameterGroup& param)
122 std::string control = param.getDefault(
"timestep.control", std::string(
"pid") );
124 const int defaultTargetIterations = 30;
125 const int defaultTargetNewtonIterations = 8;
127 const double tol = param.getDefault(
"timestep.control.tol",
double(1e-1) );
128 if( control ==
"pid" ) {
131 else if ( control ==
"pid+iteration" )
133 const int iterations = param.getDefault(
"timestep.control.targetiteration", defaultTargetIterations );
134 timeStepControl_ = TimeStepControlType(
new PIDAndIterationCountTimeStepControl( iterations, tol ) );
136 else if ( control ==
"pid+newtoniteration" )
138 const int iterations = param.getDefault(
"timestep.control.targetiteration", defaultTargetNewtonIterations );
139 timeStepControl_ = TimeStepControlType(
new PIDAndIterationCountTimeStepControl( iterations, tol ) );
142 else if ( control ==
"iterationcount" )
144 const int iterations = param.getDefault(
"timestep.control.targetiteration", defaultTargetIterations );
145 const double decayrate = param.getDefault(
"timestep.control.decayrate",
double(0.75) );
146 const double growthrate = param.getDefault(
"timestep.control.growthrate",
double(1.25) );
147 timeStepControl_ = TimeStepControlType(
new SimpleIterationCountTimeStepControl( iterations, decayrate, growthrate ) );
148 }
else if ( control ==
"hardcoded") {
149 const std::string filename = param.getDefault(
"timestep.control.filename", std::string(
"timesteps"));
150 timeStepControl_ = TimeStepControlType(
new HardcodedTimeStepControl( filename ) );
154 OPM_THROW(std::runtime_error,
"Unsupported time step control selected "<< control );
162 template <
class Solver,
class State,
class WellState>
164 step(
const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state,
const bool event )
166 return stepImpl( simulatorTimer, solver, state, well_state, event,
nullptr,
nullptr );
169 template <
class Solver,
class State,
class WellState,
class Output>
172 Solver& solver, State& state, WellState& well_state,
174 Output& outputWriter,
175 const std::vector<int>* fipnum)
177 return stepImpl( simulatorTimer, solver, state, well_state, event, &outputWriter, fipnum );
182 template <
class Solver,
class State,
class WState,
class Output >
183 SimulatorReport AdaptiveTimeStepping::
185 Solver& solver, State& state, WState& well_state,
187 Output* outputWriter,
188 const std::vector<int>* fipnum)
190 SimulatorReport report;
211 State last_state( state );
212 WState last_well_state( well_state );
221 while( ! substepTimer.done() )
224 const double dt = substepTimer.currentStepLength() ;
227 std::ostringstream ss;
228 ss <<
" Substep " << substepTimer.currentStepNum() <<
", stepsize "
229 << unit::convert::to(substepTimer.currentStepLength(), unit::day) <<
" days.";
230 OpmLog::info(ss.str());
233 SimulatorReport substepReport;
234 std::string cause_of_failure =
"";
236 substepReport = solver.step( substepTimer, state, well_state);
237 report += substepReport;
241 OpmLog::note(
"Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
244 catch (
const Opm::TooManyIterations& e) {
245 substepReport += solver.failureReport();
246 cause_of_failure =
"Solver convergence failure - Iteration limit reached";
251 catch (
const Opm::LinearSolverProblem& e) {
252 substepReport += solver.failureReport();
253 cause_of_failure =
"Linear solver convergence failure";
258 catch (
const Opm::NumericalProblem& e) {
259 substepReport += solver.failureReport();
260 cause_of_failure =
"Solver convergence failure - Numerical problem encountered";
265 catch (
const std::runtime_error& e) {
266 substepReport += solver.failureReport();
271 catch (
const Dune::ISTLError& e) {
272 substepReport += solver.failureReport();
277 catch (
const Dune::MatrixBlockError& e) {
278 substepReport += solver.failureReport();
284 if( substepReport.converged )
290 detail::SolutionTimeErrorSolverWrapper< Solver, State >
291 relativeChange( solver, last_state, state );
295 : substepReport.total_linear_iterations;
296 double dtEstimate =
timeStepControl_->computeTimeStepSize( dt, iterations, relativeChange,
297 substepTimer.simulationTimeElapsed());
300 dtEstimate = std::min( dtEstimate,
double(
max_growth_ * dt) );
311 std::ostringstream ss;
312 ss <<
" Substep summary: ";
313 if (substepReport.total_well_iterations != 0) {
314 ss <<
"well its = " << std::setw(2) << substepReport.total_well_iterations <<
", ";
316 ss <<
"newton its = " << std::setw(2) << substepReport.total_newton_iterations <<
", "
317 <<
"linearizations = " << std::setw(2) << substepReport.total_linearizations
318 <<
" (" << std::fixed << std::setprecision(3) << std::setw(6) << substepReport.assemble_time <<
" sec), "
319 <<
"linear its = " << std::setw(3) << substepReport.total_linear_iterations
320 <<
" (" << std::fixed << std::setprecision(3) << std::setw(6) << substepReport.linear_solve_time <<
" sec)";
321 OpmLog::info(ss.str());
328 if( outputWriter && !substepTimer.done() ) {
330 solver.computeFluidInPlace(state, *fipnum);
332 Opm::time::StopWatch perfTimer;
335 const auto& physicalModel = solver.model();
336 outputWriter->writeTimeStep( substepTimer, state, well_state, physicalModel, substep);
337 report.output_write_time += perfTimer.secsSinceStart();
341 substepTimer.provideTimeStepEstimate( dtEstimate );
345 last_well_state = well_state;
347 report.converged = substepTimer.done();
348 substepTimer.setLastStepFailed(
false);
353 substepTimer.setLastStepFailed(
true);
359 const auto msg = std::string(
"Solver failed to converge after cutting timestep ")
360 + std::to_string(restarts) +
" times.";
364 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
369 substepTimer.provideTimeStepEstimate( newTimeStep );
372 msg = cause_of_failure +
"\nTimestep chopped to "
373 + std::to_string(unit::convert::to( substepTimer.currentStepLength(), unit::day )) +
" days\n";
374 OpmLog::problem(msg);
378 well_state = last_well_state;
389 std::ostringstream ss;
390 substepTimer.report(ss);
391 ss <<
"Suggested next step size = " << unit::convert::to(
suggested_next_timestep_, unit::day ) <<
" (days)" << std::endl;
392 OpmLog::note(ss.str());
TimeStepControlType timeStepControl_
time step control object
Definition: AdaptiveTimeStepping.hpp:113
double suggested_next_timestep_
suggested size of next timestep
Definition: AdaptiveTimeStepping.hpp:121
const double max_time_step_
maximal allowed time step size
Definition: AdaptiveTimeStepping.hpp:117
SimulatorReport step(const SimulatorTimer &timer, Solver &solver, State &state, WellState &well_state, const bool event)
step method that acts like the solver::step method in a sub cycle of time steps
Definition: AdaptiveTimeStepping_impl.hpp:164
SimulatorReport failureReport_
statistics for the failed substeps of the last timestep
Definition: AdaptiveTimeStepping.hpp:112
double currentStepLength() const
Current step length.
Definition: SimulatorTimer.cpp:91
const double restart_factor_
factor to multiply time step with when solver fails to converge
Definition: AdaptiveTimeStepping.hpp:114
const double max_growth_
factor that limits the maximum growth of a time step
Definition: AdaptiveTimeStepping.hpp:116
const double growth_factor_
factor to multiply time step when solver recovered from failed convergence
Definition: AdaptiveTimeStepping.hpp:115
AdaptiveTimeStepping(const ParameterGroup ¶m, const bool terminal_output=true)
contructor taking parameter object
Definition: AdaptiveTimeStepping_impl.hpp:99
const int solver_restart_max_
how many restart of solver are allowed
Definition: AdaptiveTimeStepping.hpp:118
const double timestep_after_event_
suggested size of timestep after an event
Definition: AdaptiveTimeStepping.hpp:123
const bool timestep_verbose_
timestep verbosity
Definition: AdaptiveTimeStepping.hpp:120
Definition: AdaptiveTimeStepping_impl.hpp:41
RelativeChangeInterface.
Definition: TimeStepControlInterface.hpp:31
const bool solver_verbose_
solver verbosity
Definition: AdaptiveTimeStepping.hpp:119
PID controller based adaptive time step control as suggested in: Turek and Kuzmin.
Definition: TimeStepControl.hpp:74
Definition: SimulatorTimer.hpp:34
bool use_newton_iteration_
use newton iteration count for adaptive time step control
Definition: AdaptiveTimeStepping.hpp:124
double relativeChange() const
return || u^n+1 - u^n || / || u^n+1 ||
Definition: AdaptiveTimeStepping_impl.hpp:56
bool full_timestep_initially_
beginning with the size of the time step from data file
Definition: AdaptiveTimeStepping.hpp:122