DUNE-FEM (unstable)

timestepcontrol.hh
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
10
11//- dune-fem includes
12#include <dune/fem/io/parameter.hh>
13#include <dune/fem/solver/timeprovider.hh>
14
15namespace DuneODE
16{
17
18 // ImplicitRungeKuttaSolverParameters
19 // ----------------------------------
20
21 struct ImplicitRungeKuttaSolverParameters
22 : public Dune::Fem::LocalParameter< ImplicitRungeKuttaSolverParameters, ImplicitRungeKuttaSolverParameters >
23 {
24 enum { noVerbosity = 0, noConvergenceVerbosity = 1,
25 cflVerbosity = 2, fullVerbosity = 3 };
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
40 Dune::Fem::ParameterReader parameter_;
41
42 public:
43 ImplicitRungeKuttaSolverParameters ( const std::string keyPrefix, const Dune::Fem::ParameterReader &parameter = Dune::Fem::Parameter::container() )
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 ) ),
48 parameter_( parameter )
49 {}
50
51 ImplicitRungeKuttaSolverParameters ( const Dune::Fem::ParameterReader &parameter = Dune::Fem::Parameter::container() )
52 : ImplicitRungeKuttaSolverParameters( "fem.ode.", parameter )
53 {}
54
55 // destructor (virtual)
56 virtual ~ImplicitRungeKuttaSolverParameters() {}
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
165 class ImplicitRungeKuttaTimeStepControl
166 {
167 typedef ImplicitRungeKuttaTimeStepControl ThisType;
168
169 public:
170 typedef Dune::Fem::TimeProviderBase TimeProviderType;
171 typedef ImplicitRungeKuttaSolverParameters ParameterType;
172
173 ImplicitRungeKuttaTimeStepControl ( TimeProviderType &timeProvider, const ParameterType &parameters )
174 : timeProvider_( timeProvider ),
175 parameters_( parameters.clone() ),
176 cfl_( parameters_->cflStart() ),
177 cflMax_( parameters_->cflMax() ),
178 verbose_( parameters_->verbose() ),
179 initialized_( false )
180 {}
181
182 explicit ImplicitRungeKuttaTimeStepControl ( TimeProviderType &timeProvider, const Dune::Fem::ParameterReader &parameter = Dune::Fem::Parameter::container() )
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_ );
198 timeProvider_.provideTimeStepEstimate( estimate );
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
216 if( (verbose_ >= ImplicitRungeKuttaSolverParameters::noConvergenceVerbosity) && (Dune::Fem::MPIManager::rank() == 0) )
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
225 timeProvider_.provideTimeStepEstimate( factor * timeStepSize() );
226 timeProvider_.invalidateTimeStep();
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
247 if( (cfl_ != oldCfl) && (verbose_ >= ImplicitRungeKuttaSolverParameters::cflVerbosity) && (Dune::Fem::MPIManager::rank() == 0) )
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
265 TimeProviderType &timeProvider_;
266 std::shared_ptr< const ParameterType > parameters_;
267 double cfl_, cflMax_;
268 int verbose_;
269 bool initialized_;
270 };
271
272
273
274 // PIDTimeStepControl
275 // ------------------
276
287 class PIDTimeStepControl : public ImplicitRungeKuttaTimeStepControl
288 {
289 typedef PIDTimeStepControl ThisType;
290 typedef ImplicitRungeKuttaTimeStepControl BaseType;
291
292 protected:
293 using BaseType :: initialized_;
294 using BaseType :: cfl_;
295 using BaseType :: parameters ;
296 public:
298 typedef typename BaseType::ParameterType ParameterType;
299
300 PIDTimeStepControl ( TimeProviderType &timeProvider, const ParameterType &parameters )
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
312 PIDTimeStepControl ( TimeProviderType &timeProvider, const Dune::Fem::ParameterReader &parameter = Dune::Fem::Parameter::container() )
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;
345 timeProvider_.provideTimeStepEstimate( dtEst );
346
347 if( (verbose_ >= ImplicitRungeKuttaSolverParameters::cflVerbosity) && (Dune::Fem::MPIManager::rank() == 0) )
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
PID time step control.
Definition: timestepcontrol.hh:288
general base for time providers
Definition: timeprovider.hh:36
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
DErrType derr(std::cerr)
Stream for error messages.
Definition: stdstreams.hh:196
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)