dune-fem  2.8-git
timestepcontrol.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_TIMESTEPCONTROL_HH
2 #define DUNE_FEM_SOLVER_RUNGEKUTTA_TIMESTEPCONTROL_HH
3 
4 //- system includes
5 #include <cassert>
6 #include <memory>
7 
8 //- dune-common includes
9 #include <dune/common/exceptions.hh>
10 
11 //- dune-fem includes
12 #include <dune/fem/io/parameter.hh>
14 
15 namespace DuneODE
16 {
17 
18  // ImplicitRungeKuttaSolverParameters
19  // ----------------------------------
20 
22  : public Dune::Fem::LocalParameter< ImplicitRungeKuttaSolverParameters, ImplicitRungeKuttaSolverParameters >
23  {
26 
27  protected:
28  // key prefix, default is fem.ode (can be overloaded by user)
29  const std::string keyPrefix_;
30 
31  // number of minimal iterations that the linear solver should do
32  // if the number of iterations done is smaller then the cfl number is increased
33  const int minIter_;
34  // number of maximal iterations that the linear solver should do
35  // if the number of iterations larger then the cfl number is decreased
36  const int maxIter_;
37  // factor for cfl number on increase (decrease is 0.5)
38  const double sigma_;
39 
41 
42  public:
44  : keyPrefix_( keyPrefix ),
45  minIter_( parameter.getValue< int >( keyPrefix_ + "miniterations", 14 ) ),
46  maxIter_( parameter.getValue< int >( keyPrefix_ + "maxiterations" , 16 ) ),
47  sigma_( parameter.getValue< double >( keyPrefix_ + "cflincrease" , 1.1 ) ),
49  {}
50 
53  {}
54 
55  // destructor (virtual)
57 
58  const Dune::Fem::ParameterReader &parameter () const { return parameter_; }
59 
62  virtual double tolerance () const
63  {
64  return parameter().getValue< double >( keyPrefix_ + "tolerance" , 1e-6 );
65  }
66 
67  virtual int iterations() const
68  {
69  return parameter().getValue< int >( keyPrefix_ + "iterations" , 1000 );
70  }
71 
73  virtual int verbose () const
74  {
75  static const std::string verboseTypeTable[] = { "none", "noconv", "cfl", "full" };
76  return parameter().getEnum( keyPrefix_ + "verbose", verboseTypeTable, 0 );
77  }
78 
79  virtual double cflStart () const
80  {
81  return parameter().getValue< double >( keyPrefix_ + "cflStart", 1 );
82  }
83 
84  virtual double cflMax () const
85  {
86  return parameter().getValue< double >( keyPrefix_ + "cflMax" , std::numeric_limits< double >::max() );
87  }
88 
89  double initialDeltaT ( double dt ) const
90  {
91  return std::min( parameter().getValue< double >( keyPrefix_ + "initialdt", 987654321 ), dt );
92  }
93 
104  virtual bool cflFactor( const double imOpTimeStepEstimate,
105  const double exOpTimeStepEstimate,
106  const int numberOfLinearIterations,
107  bool converged,
108  double &factor) const
109  {
110  const int iter = numberOfLinearIterations;
111  factor = 1.;
112  bool changed = false;
113  if (converged)
114  {
115  if (iter < minIter_)
116  {
117  if( imOpTimeStepEstimate <= exOpTimeStepEstimate )
118  {
119  factor = sigma_;
120  changed = true;
121  }
122  }
123  else if (iter > maxIter_)
124  {
125  factor = (double)maxIter_/(sigma_*(double)iter);
126  changed = true;
127  }
128  }
129  else
130  {
131  factor = 0.5;
132  changed = true;
133  }
134  return changed;
135  }
136 
137  virtual void initTimeStepEstimate ( const double dtEstExpl, const double dtEstImpl, double &dtEst, double &cfl ) const
138  {
139  // initial time step already set to explicit time step
140  dtEst = dtEstExpl;
141 
142  // heuristics for initial CFL number
143  cfl = 1.0;
144  if( (dtEstImpl > 0) && (dtEstExpl > dtEstImpl) )
145  cfl = dtEstExpl / dtEstImpl;
146  }
147 
148  // return number of max linear iterations per newton step
149  virtual int maxLinearIterations () const { return maxIter_; }
150 
152  virtual int selectedSolver( const int order ) const
153  {
154  const std::string names [] = { "ImplicitEuler", "CrankNicolson", "DIRK23", "DIRK34", "SDIRK22" };
155  // by default select according to order
156  return parameter().getEnum( keyPrefix_ + "solvername", names, order-1 ) + 1;
157  }
158  };
159 
160 
161 
162  // ImplicitRungeKuttaTimeStepControl
163  // ---------------------------------
164 
166  {
168 
169  public:
172 
174  : timeProvider_( timeProvider ),
175  parameters_( parameters.clone() ),
176  cfl_( parameters_->cflStart() ),
177  cflMax_( parameters_->cflMax() ),
178  verbose_( parameters_->verbose() ),
179  initialized_( false )
180  {}
181 
183  : timeProvider_( timeProvider ),
184  parameters_( std::make_shared<ParameterType>( parameter ) ),
185  cfl_( parameters_->cflStart() ),
186  cflMax_( parameters_->cflMax() ),
187  verbose_( parameters_->verbose() ),
188  initialized_( false )
189  {}
190 
191  double time () const { return timeProvider_.time(); }
192  double timeStepSize () const { return timeProvider_.deltaT(); }
193 
194  void initialTimeStepSize ( double helmholtzEstimate, double sourceTermEstimate )
195  {
196  double estimate = std::numeric_limits< double >::max();
197  parameters().initTimeStepEstimate( sourceTermEstimate, helmholtzEstimate, estimate, cfl_ );
199  initialized_ = true;
200  }
201 
202  template< class Monitor >
203  void reduceTimeStep ( double helmholtzEstimate, double sourceTermEstimate, const Monitor &monitor )
204  {
205  if( !initialized_ )
206  DUNE_THROW( Dune::InvalidStateException, "ImplicitRungeKuttaSolver must be initialized before first solve." );
207 
208  double factor( 1 );
209  parameters().cflFactor( helmholtzEstimate, sourceTermEstimate, monitor.linearSolverIterations_, false, factor );
210 
211  if( !((factor >= std::numeric_limits< double >::min()) && (factor < 1.0)) )
212  DUNE_THROW( Dune::InvalidStateException, "invalid cfl factor: " << factor );
213 
214  cfl_ *= factor;
215 
217  {
218  Dune::derr << "Implicit Runge-Kutta step failed to converge." << std::endl;
219  Dune::derr << "New cfl number: " << cfl_
220  << ", iterations per time step: ILS = " << monitor.linearSolverIterations_
221  << ", INLS = " << monitor.newtonIterations_
222  << std::endl;
223  }
224 
227  }
228 
229  template< class Monitor >
230  void timeStepEstimate ( double helmholtzEstimate, double sourceTermEstimate, const Monitor &monitor )
231  {
232  if( !initialized_ )
233  DUNE_THROW( Dune::InvalidStateException, "ImplicitRungeKuttaSolver must be initialized before first solve." );
234 
235  double factor( 1 );
236  // true means converged, which is always true since this function is only called
237  // when the implicit solver did converge
238  parameters().cflFactor( helmholtzEstimate, sourceTermEstimate, monitor.linearSolverIterations_, true, factor );
239  if( !((factor >= std::numeric_limits< double >::min()) && (factor <= std::numeric_limits< double >::max())) )
240  DUNE_THROW( Dune::InvalidStateException, "invalid cfl factor: " << factor );
241 
242  const double oldCfl = cfl_;
243  cfl_ = std::min( cflMax_, factor * cfl_ );
244 
245  timeProvider_.provideTimeStepEstimate( std::min( sourceTermEstimate, cfl_ * helmholtzEstimate ) );
246 
248  {
249  Dune::derr << "New cfl number: " << cfl_
250  << ", iterations per time step: ILS = " << monitor.linearSolverIterations_
251  << ", INLS = " << monitor.newtonIterations_
252  << std::endl;
253  }
254  }
255 
256  bool computeError () const { return false; }
257 
258  protected:
259  const ParameterType &parameters () const
260  {
261  assert( parameters_ );
262  return *parameters_;
263  }
264 
266  std::shared_ptr< const ParameterType > parameters_;
267  double cfl_, cflMax_;
268  int verbose_;
270  };
271 
272 
273 
274  // PIDTimeStepControl
275  // ------------------
276 
288  {
291 
292  protected:
294  using BaseType :: cfl_;
295  using BaseType :: parameters ;
296  public:
299 
301  : BaseType( timeProvider, parameters ),
302  errors_(),
303  tol_( 1e-3 )
304  {
305  if( parameters.parameter().getValue("fem.ode.pidcontrol", bool(false) ) )
306  {
307  tol_ = parameters.parameter().getValue("fem.ode.pidtolerance", tol_ );
308  errors_.resize( 3, tol_ );
309  }
310  }
311 
313  : BaseType( timeProvider, parameter ),
314  errors_(),
315  tol_( 1e-3 )
316  {
317  if( parameter.getValue("fem.ode.pidcontrol", bool(false) ) )
318  {
319  tol_ = parameter.getValue("fem.ode.pidtolerance", tol_ );
320  errors_.resize( 3, tol_ );
321  }
322  }
323 
324  bool computeError () const { return ! errors_.empty() ; }
325 
326  template< class Monitor >
327  void timeStepEstimate ( double helmholtzEstimate, double sourceTermEstimate, const Monitor &monitor )
328  {
329  if( !initialized_ )
330  DUNE_THROW( Dune::InvalidStateException, "ImplicitRungeKuttaSolver must be initialized before first solve." );
331 
332  if( computeError() ) // use pid control
333  {
334  cfl_ = 1.0; // reset cfl for next reduceTimeStep
335  double dtEst = pidTimeStepControl( std::min( sourceTermEstimate, helmholtzEstimate ), monitor );
336  /*
337  const int targetIterations = parameters().maxLinearIterations();
338  if( monitor.linearSolverIterations_ > targetIterations &&
339  targetIterations > 0 )
340  {
341  dtEst *= double( targetIterations ) / double(monitor.linearSolverIterations_);
342  }
343  */
344  std::cout << "Set dt = " << dtEst << std::endl;
346 
348  {
349  Dune::derr << "New dt: " << dtEst
350  << ", iterations per time step: ILS = " << monitor.linearSolverIterations_
351  << ", INLS = " << monitor.newtonIterations_
352  << std::endl;
353  }
354  }
355  else
356  {
357  BaseType::timeStepEstimate( helmholtzEstimate, sourceTermEstimate, monitor );
358  }
359  }
360 
361  template < class Monitor >
362  double pidTimeStepControl( const double dt, const Monitor& monitor )
363  {
364  // get error || u^n - u^n+1 || / || u^n+1 || from monitor
365  const double error = monitor.error_;
366  std::cout << error << " error " << std::endl;
367  if( std::abs( error ) < 1e-12 ) return 10. * dt;
368 
369  // shift errors
370  for( int i=0; i<2; ++i )
371  {
372  errors_[ i ] = errors_[i+1];
373  }
374 
375  // store new error
376  errors_[ 2 ] = error ;
377 
378  if( error > tol_ )
379  {
380  // adjust dt by given tolerance
381  const double newDt = dt * tol_ / error;
382  return newDt;
383  }
384  else if( error > 1e-12 )
385  {
386  // values taking from turek time stepping paper
387  const double kP = 0.075 ;
388  const double kI = 0.175 ;
389  const double kD = 0.01 ;
390  const double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) *
391  std::pow( tol_ / errors_[ 2 ], kI ) *
392  std::pow( errors_[0]*errors_[0]/errors_[ 1 ]/errors_[ 2 ], kD ));
393  std::cout << "newDt = " << newDt << std::endl;
394  return newDt;
395  }
396 
397  return dt ;
398  }
399 
400  protected:
401  std::vector< double > errors_;
402  double tol_;
403  };
404 
405 } // namespace DuneODE
406 
407 #endif // #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_TIMESTEPCONTROL_HH
double pow(const Double &a, const double b)
Definition: double.hh:881
Double abs(const Double &a)
Definition: double.hh:871
static constexpr T max(T a)
Definition: utility.hh:77
static constexpr T min(T a)
Definition: utility.hh:93
Definition: multistep.hh:17
int getEnum(const std::string &key, const std::string(&values)[n]) const
Definition: reader.hh:225
T getValue(const std::string &key) const
get mandatory parameter
Definition: reader.hh:159
static ParameterContainer & container()
Definition: io/parameter.hh:193
Definition: io/parameter.hh:551
static int rank()
Definition: mpimanager.hh:155
Definition: timestepcontrol.hh:23
virtual double cflMax() const
Definition: timestepcontrol.hh:84
double initialDeltaT(double dt) const
Definition: timestepcontrol.hh:89
virtual int maxLinearIterations() const
Definition: timestepcontrol.hh:149
virtual double tolerance() const
tolerance for the non-linear solver (should be larger than the tolerance for the linear solver
Definition: timestepcontrol.hh:62
ImplicitRungeKuttaSolverParameters(const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: timestepcontrol.hh:51
virtual int verbose() const
verbosity level ( none, noconv, cfl, full )
Definition: timestepcontrol.hh:73
virtual ~ImplicitRungeKuttaSolverParameters()
Definition: timestepcontrol.hh:56
const int maxIter_
Definition: timestepcontrol.hh:36
virtual double cflStart() const
Definition: timestepcontrol.hh:79
ImplicitRungeKuttaSolverParameters(const std::string keyPrefix, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: timestepcontrol.hh:43
virtual void initTimeStepEstimate(const double dtEstExpl, const double dtEstImpl, double &dtEst, double &cfl) const
Definition: timestepcontrol.hh:137
const int minIter_
Definition: timestepcontrol.hh:33
@ cflVerbosity
Definition: timestepcontrol.hh:25
@ noConvergenceVerbosity
Definition: timestepcontrol.hh:24
@ fullVerbosity
Definition: timestepcontrol.hh:25
@ noVerbosity
Definition: timestepcontrol.hh:24
const std::string keyPrefix_
Definition: timestepcontrol.hh:29
const double sigma_
Definition: timestepcontrol.hh:38
virtual bool cflFactor(const double imOpTimeStepEstimate, const double exOpTimeStepEstimate, const int numberOfLinearIterations, bool converged, double &factor) const
return multiplication factor for the current cfl number
Definition: timestepcontrol.hh:104
virtual int iterations() const
Definition: timestepcontrol.hh:67
const Dune::Fem::ParameterReader & parameter() const
Definition: timestepcontrol.hh:58
Dune::Fem::ParameterReader parameter_
Definition: timestepcontrol.hh:40
virtual int selectedSolver(const int order) const
return number of selected solver (default = order of solver)
Definition: timestepcontrol.hh:152
Definition: timestepcontrol.hh:166
TimeProviderType & timeProvider_
Definition: timestepcontrol.hh:265
double cfl_
Definition: timestepcontrol.hh:267
ImplicitRungeKuttaSolverParameters ParameterType
Definition: timestepcontrol.hh:171
void initialTimeStepSize(double helmholtzEstimate, double sourceTermEstimate)
Definition: timestepcontrol.hh:194
ImplicitRungeKuttaTimeStepControl(TimeProviderType &timeProvider, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: timestepcontrol.hh:182
void reduceTimeStep(double helmholtzEstimate, double sourceTermEstimate, const Monitor &monitor)
Definition: timestepcontrol.hh:203
const ParameterType & parameters() const
Definition: timestepcontrol.hh:259
Dune::Fem::TimeProviderBase TimeProviderType
Definition: timestepcontrol.hh:170
bool computeError() const
Definition: timestepcontrol.hh:256
void timeStepEstimate(double helmholtzEstimate, double sourceTermEstimate, const Monitor &monitor)
Definition: timestepcontrol.hh:230
std::shared_ptr< const ParameterType > parameters_
Definition: timestepcontrol.hh:266
double timeStepSize() const
Definition: timestepcontrol.hh:192
int verbose_
Definition: timestepcontrol.hh:268
double time() const
Definition: timestepcontrol.hh:191
double cflMax_
Definition: timestepcontrol.hh:267
ImplicitRungeKuttaTimeStepControl(TimeProviderType &timeProvider, const ParameterType &parameters)
Definition: timestepcontrol.hh:173
bool initialized_
Definition: timestepcontrol.hh:269
PID time step control.
Definition: timestepcontrol.hh:288
double tol_
Definition: timestepcontrol.hh:402
PIDTimeStepControl(TimeProviderType &timeProvider, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: timestepcontrol.hh:312
PIDTimeStepControl(TimeProviderType &timeProvider, const ParameterType &parameters)
Definition: timestepcontrol.hh:300
BaseType::ParameterType ParameterType
Definition: timestepcontrol.hh:298
Dune::Fem::TimeProviderBase TimeProviderType
Definition: timestepcontrol.hh:297
double pidTimeStepControl(const double dt, const Monitor &monitor)
Definition: timestepcontrol.hh:362
std::vector< double > errors_
Definition: timestepcontrol.hh:401
void timeStepEstimate(double helmholtzEstimate, double sourceTermEstimate, const Monitor &monitor)
Definition: timestepcontrol.hh:327
bool computeError() const
Definition: timestepcontrol.hh:324
general base for time providers
Definition: timeprovider.hh:36
double time() const
obtain the current time
Definition: timeprovider.hh:94
void provideTimeStepEstimate(const double dtEstimate)
set time step estimate to minimum of given value and internal time step estiamte
Definition: timeprovider.hh:142
double deltaT() const
obtain the size of the current time step
Definition: timeprovider.hh:113
void invalidateTimeStep()
count current time step a not valid
Definition: timeprovider.hh:158