00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
00020 #define OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
00021
00022 #include <iostream>
00023 #include <string>
00024 #include <utility>
00025
00026 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
00027 #include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
00028 #include <opm/simulators/timestepping/TimeStepControl.hpp>
00029 #include <opm/core/utility/StopWatch.hpp>
00030 #include <opm/common/Exceptions.hpp>
00031 #include <opm/common/OpmLog/OpmLog.hpp>
00032 #include <dune/istl/istlexception.hh>
00033 #include <dune/istl/ilu.hh>
00034 #include <opm/parser/eclipse/EclipseState/Schedule/Tuning.hpp>
00035
00036 namespace Opm {
00037
00038 namespace detail
00039 {
00040 template <class Solver, class State>
00041 class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface
00042 {
00043 const Solver& solver_;
00044 const State& previous_;
00045 const State& current_;
00046 public:
00047 SolutionTimeErrorSolverWrapper( const Solver& solver,
00048 const State& previous,
00049 const State& current )
00050 : solver_( solver ),
00051 previous_( previous ),
00052 current_( current )
00053 {}
00054
00056 double relativeChange() const
00057 {
00058 return solver_.model().relativeChange( previous_, current_ );
00059 }
00060 };
00061
00062 template<class E>
00063 void logException(const E& exception, bool verbose)
00064 {
00065 if( verbose )
00066 {
00067 std::string message;
00068 message = "Caught Exception: ";
00069 message += exception.what();
00070 OpmLog::debug(message);
00071 }
00072 }
00073 }
00074
00075
00076
00077 inline AdaptiveTimeStepping::AdaptiveTimeStepping( const Tuning& tuning,
00078 size_t time_step,
00079 const ParameterGroup& param,
00080 const bool terminal_output )
00081 : timeStepControl_()
00082 , restart_factor_( tuning.getTSFCNV(time_step) )
00083 , growth_factor_(tuning.getTFDIFF(time_step) )
00084 , max_growth_( tuning.getTSFMAX(time_step) )
00085
00086 , max_time_step_( tuning.getTSMAXZ(time_step) )
00087 , solver_restart_max_( param.getDefault("solver.restart", int(10) ) )
00088 , solver_verbose_( param.getDefault("solver.verbose", bool(true) ) && terminal_output )
00089 , timestep_verbose_( param.getDefault("timestep.verbose", bool(true) ) && terminal_output )
00090 , suggested_next_timestep_( tuning.getTSINIT(time_step) )
00091 , full_timestep_initially_( param.getDefault("full_timestep_initially", bool(false) ) )
00092 , timestep_after_event_( tuning.getTMAXWC(time_step))
00093 , use_newton_iteration_(false)
00094 {
00095 init(param);
00096
00097 }
00098
00099 inline AdaptiveTimeStepping::AdaptiveTimeStepping( const ParameterGroup& param,
00100 const bool terminal_output )
00101 : timeStepControl_()
00102 , restart_factor_( param.getDefault("solver.restartfactor", double(0.33) ) )
00103 , growth_factor_( param.getDefault("solver.growthfactor", double(2) ) )
00104 , max_growth_( param.getDefault("timestep.control.maxgrowth", double(3.0) ) )
00105
00106 , max_time_step_( unit::convert::from(param.getDefault("timestep.max_timestep_in_days", 365.0 ), unit::day) )
00107 , solver_restart_max_( param.getDefault("solver.restart", int(10) ) )
00108 , solver_verbose_( param.getDefault("solver.verbose", bool(true) ) && terminal_output )
00109 , timestep_verbose_( param.getDefault("timestep.verbose", bool(true) ) && terminal_output )
00110 , suggested_next_timestep_( unit::convert::from(param.getDefault("timestep.initial_timestep_in_days", -1.0 ), unit::day) )
00111 , full_timestep_initially_( param.getDefault("full_timestep_initially", bool(false) ) )
00112 , timestep_after_event_( unit::convert::from(param.getDefault("timestep.timestep_in_days_after_event", -1.0 ), unit::day))
00113 , use_newton_iteration_(false)
00114 {
00115 init(param);
00116 }
00117
00118 inline void AdaptiveTimeStepping::
00119 init(const ParameterGroup& param)
00120 {
00121
00122 std::string control = param.getDefault("timestep.control", std::string("pid") );
00123
00124 const int defaultTargetIterations = 30;
00125 const int defaultTargetNewtonIterations = 8;
00126
00127 const double tol = param.getDefault("timestep.control.tol", double(1e-1) );
00128 if( control == "pid" ) {
00129 timeStepControl_ = TimeStepControlType( new PIDTimeStepControl( tol ) );
00130 }
00131 else if ( control == "pid+iteration" )
00132 {
00133 const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
00134 timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
00135 }
00136 else if ( control == "pid+newtoniteration" )
00137 {
00138 const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetNewtonIterations );
00139 timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
00140 use_newton_iteration_ = true;
00141 }
00142 else if ( control == "iterationcount" )
00143 {
00144 const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
00145 const double decayrate = param.getDefault("timestep.control.decayrate", double(0.75) );
00146 const double growthrate = param.getDefault("timestep.control.growthrate", double(1.25) );
00147 timeStepControl_ = TimeStepControlType( new SimpleIterationCountTimeStepControl( iterations, decayrate, growthrate ) );
00148 } else if ( control == "hardcoded") {
00149 const std::string filename = param.getDefault("timestep.control.filename", std::string("timesteps"));
00150 timeStepControl_ = TimeStepControlType( new HardcodedTimeStepControl( filename ) );
00151
00152 }
00153 else
00154 OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control );
00155
00156
00157 assert( growth_factor_ >= 1.0 );
00158 }
00159
00160
00161
00162 template <class Solver, class State, class WellState>
00163 SimulatorReport AdaptiveTimeStepping::
00164 step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state, const bool event )
00165 {
00166 return stepImpl( simulatorTimer, solver, state, well_state, event, nullptr, nullptr );
00167 }
00168
00169 template <class Solver, class State, class WellState, class Output>
00170 SimulatorReport AdaptiveTimeStepping::
00171 step( const SimulatorTimer& simulatorTimer,
00172 Solver& solver, State& state, WellState& well_state,
00173 const bool event,
00174 Output& outputWriter,
00175 const std::vector<int>* fipnum)
00176 {
00177 return stepImpl( simulatorTimer, solver, state, well_state, event, &outputWriter, fipnum );
00178 }
00179
00180
00181
00182 template <class Solver, class State, class WState, class Output >
00183 SimulatorReport AdaptiveTimeStepping::
00184 stepImpl( const SimulatorTimer& simulatorTimer,
00185 Solver& solver, State& state, WState& well_state,
00186 const bool event,
00187 Output* outputWriter,
00188 const std::vector<int>* fipnum)
00189 {
00190 SimulatorReport report;
00191 const double timestep = simulatorTimer.currentStepLength();
00192
00193
00194 if( suggested_next_timestep_ < 0 ) {
00195 suggested_next_timestep_ = restart_factor_ * timestep;
00196 }
00197
00198 if (full_timestep_initially_) {
00199 suggested_next_timestep_ = timestep;
00200 }
00201
00202
00203 if (event && timestep_after_event_ > 0) {
00204 suggested_next_timestep_ = timestep_after_event_;
00205 }
00206
00207
00208 AdaptiveSimulatorTimer substepTimer( simulatorTimer, suggested_next_timestep_, max_time_step_ );
00209
00210
00211 State last_state( state );
00212 WState last_well_state( well_state );
00213
00214
00215 failureReport_ = SimulatorReport();
00216
00217
00218 int restarts = 0;
00219
00220
00221 while( ! substepTimer.done() )
00222 {
00223
00224 const double dt = substepTimer.currentStepLength() ;
00225 if( timestep_verbose_ )
00226 {
00227 std::ostringstream ss;
00228 ss <<" Substep " << substepTimer.currentStepNum() << ", stepsize "
00229 << unit::convert::to(substepTimer.currentStepLength(), unit::day) << " days.";
00230 OpmLog::info(ss.str());
00231 }
00232
00233 SimulatorReport substepReport;
00234 std::string cause_of_failure = "";
00235 try {
00236 substepReport = solver.step( substepTimer, state, well_state);
00237 report += substepReport;
00238
00239 if( solver_verbose_ ) {
00240
00241 OpmLog::note("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
00242 }
00243 }
00244 catch (const Opm::TooManyIterations& e) {
00245 substepReport += solver.failureReport();
00246 cause_of_failure = "Solver convergence failure - Iteration limit reached";
00247
00248 detail::logException(e, solver_verbose_);
00249
00250 }
00251 catch (const Opm::LinearSolverProblem& e) {
00252 substepReport += solver.failureReport();
00253 cause_of_failure = "Linear solver convergence failure";
00254
00255 detail::logException(e, solver_verbose_);
00256
00257 }
00258 catch (const Opm::NumericalProblem& e) {
00259 substepReport += solver.failureReport();
00260 cause_of_failure = "Solver convergence failure - Numerical problem encountered";
00261
00262 detail::logException(e, solver_verbose_);
00263
00264 }
00265 catch (const std::runtime_error& e) {
00266 substepReport += solver.failureReport();
00267
00268 detail::logException(e, solver_verbose_);
00269
00270 }
00271 catch (const Dune::ISTLError& e) {
00272 substepReport += solver.failureReport();
00273
00274 detail::logException(e, solver_verbose_);
00275
00276 }
00277 catch (const Dune::MatrixBlockError& e) {
00278 substepReport += solver.failureReport();
00279
00280 detail::logException(e, solver_verbose_);
00281
00282 }
00283
00284 if( substepReport.converged )
00285 {
00286
00287 ++substepTimer;
00288
00289
00290 detail::SolutionTimeErrorSolverWrapper< Solver, State >
00291 relativeChange( solver, last_state, state );
00292
00293
00294 const int iterations = use_newton_iteration_ ? substepReport.total_newton_iterations
00295 : substepReport.total_linear_iterations;
00296 double dtEstimate = timeStepControl_->computeTimeStepSize( dt, iterations, relativeChange,
00297 substepTimer.simulationTimeElapsed());
00298
00299
00300 dtEstimate = std::min( dtEstimate, double(max_growth_ * dt) );
00301
00302
00303 if( restarts > 0 ) {
00304 dtEstimate = std::min( growth_factor_ * dt, dtEstimate );
00305
00306 restarts = 0;
00307 }
00308
00309 if( timestep_verbose_ )
00310 {
00311 std::ostringstream ss;
00312 ss << " Substep summary: ";
00313 if (substepReport.total_well_iterations != 0) {
00314 ss << "well its = " << std::setw(2) << substepReport.total_well_iterations << ", ";
00315 }
00316 ss << "newton its = " << std::setw(2) << substepReport.total_newton_iterations << ", "
00317 << "linearizations = " << std::setw(2) << substepReport.total_linearizations
00318 << " (" << std::fixed << std::setprecision(3) << std::setw(6) << substepReport.assemble_time << " sec), "
00319 << "linear its = " << std::setw(3) << substepReport.total_linear_iterations
00320 << " (" << std::fixed << std::setprecision(3) << std::setw(6) << substepReport.linear_solve_time << " sec)";
00321 OpmLog::info(ss.str());
00322 }
00323
00324
00325
00326
00327
00328 if( outputWriter && !substepTimer.done() ) {
00329 if (fipnum) {
00330 solver.computeFluidInPlace(state, *fipnum);
00331 }
00332 Opm::time::StopWatch perfTimer;
00333 perfTimer.start();
00334 bool substep = true;
00335 const auto& physicalModel = solver.model();
00336 outputWriter->writeTimeStep( substepTimer, state, well_state, physicalModel, substep);
00337 report.output_write_time += perfTimer.secsSinceStart();
00338 }
00339
00340
00341 substepTimer.provideTimeStepEstimate( dtEstimate );
00342
00343
00344 last_state = state ;
00345 last_well_state = well_state;
00346
00347 report.converged = substepTimer.done();
00348 substepTimer.setLastStepFailed(false);
00349
00350 }
00351 else
00352 {
00353 substepTimer.setLastStepFailed(true);
00354
00355 failureReport_ += substepReport;
00356
00357
00358 if( restarts >= solver_restart_max_ ) {
00359 const auto msg = std::string("Solver failed to converge after cutting timestep ")
00360 + std::to_string(restarts) + " times.";
00361 if (solver_verbose_) {
00362 OpmLog::error(msg);
00363 }
00364 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
00365 }
00366
00367 const double newTimeStep = restart_factor_ * dt;
00368
00369 substepTimer.provideTimeStepEstimate( newTimeStep );
00370 if( solver_verbose_ ) {
00371 std::string msg;
00372 msg = cause_of_failure + "\nTimestep chopped to "
00373 + std::to_string(unit::convert::to( substepTimer.currentStepLength(), unit::day )) + " days\n";
00374 OpmLog::problem(msg);
00375 }
00376
00377 state = last_state;
00378 well_state = last_well_state;
00379
00380 ++restarts;
00381 }
00382 }
00383
00384
00385
00386 suggested_next_timestep_ = substepTimer.currentStepLength();
00387 if( timestep_verbose_ )
00388 {
00389 std::ostringstream ss;
00390 substepTimer.report(ss);
00391 ss << "Suggested next step size = " << unit::convert::to( suggested_next_timestep_, unit::day ) << " (days)" << std::endl;
00392 OpmLog::note(ss.str());
00393 }
00394
00395 if( ! std::isfinite( suggested_next_timestep_ ) ) {
00396 suggested_next_timestep_ = timestep;
00397 }
00398 return report;
00399 }
00400 }
00401
00402 #endif