1#ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_TIMESTEPCONTROL_HH
2#define DUNE_FEM_SOLVER_RUNGEKUTTA_TIMESTEPCONTROL_HH
12#include <dune/fem/io/parameter.hh>
13#include <dune/fem/solver/timeprovider.hh>
21 struct ImplicitRungeKuttaSolverParameters
22 :
public Dune::Fem::LocalParameter< ImplicitRungeKuttaSolverParameters, ImplicitRungeKuttaSolverParameters >
24 enum { noVerbosity = 0, noConvergenceVerbosity = 1,
25 cflVerbosity = 2, fullVerbosity = 3 };
29 const std::string keyPrefix_;
40 Dune::Fem::ParameterReader parameter_;
43 ImplicitRungeKuttaSolverParameters (
const std::string keyPrefix,
const Dune::Fem::ParameterReader ¶meter = 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 )
51 ImplicitRungeKuttaSolverParameters (
const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
52 : ImplicitRungeKuttaSolverParameters(
"fem.ode.", parameter )
56 virtual ~ImplicitRungeKuttaSolverParameters() {}
58 const Dune::Fem::ParameterReader ¶meter ()
const {
return parameter_; }
62 virtual double tolerance ()
const
64 return parameter().getValue<
double >( keyPrefix_ +
"tolerance" , 1e-6 );
67 virtual int iterations()
const
69 return parameter().getValue<
int >( keyPrefix_ +
"iterations" , 1000 );
73 virtual int verbose ()
const
75 static const std::string verboseTypeTable[] = {
"none",
"noconv",
"cfl",
"full" };
76 return parameter().getEnum( keyPrefix_ +
"verbose", verboseTypeTable, 0 );
79 virtual double cflStart ()
const
81 return parameter().getValue<
double >( keyPrefix_ +
"cflStart", 1 );
84 virtual double cflMax ()
const
89 double initialDeltaT (
double dt )
const
91 return std::min( parameter().getValue< double >( keyPrefix_ +
"initialdt", 987654321 ), dt );
104 virtual bool cflFactor(
const double imOpTimeStepEstimate,
105 const double exOpTimeStepEstimate,
106 const int numberOfLinearIterations,
108 double &factor)
const
110 const int iter = numberOfLinearIterations;
112 bool changed =
false;
117 if( imOpTimeStepEstimate <= exOpTimeStepEstimate )
123 else if (iter > maxIter_)
125 factor = (double)maxIter_/(sigma_*(
double)iter);
137 virtual void initTimeStepEstimate (
const double dtEstExpl,
const double dtEstImpl,
double &dtEst,
double &cfl )
const
144 if( (dtEstImpl > 0) && (dtEstExpl > dtEstImpl) )
145 cfl = dtEstExpl / dtEstImpl;
149 virtual int maxLinearIterations ()
const {
return maxIter_; }
152 virtual int selectedSolver(
const int order )
const
154 const std::string names [] = {
"ImplicitEuler",
"CrankNicolson",
"DIRK23",
"DIRK34",
"SDIRK22" };
156 return parameter().getEnum( keyPrefix_ +
"solvername", names, order-1 ) + 1;
165 class ImplicitRungeKuttaTimeStepControl
167 typedef ImplicitRungeKuttaTimeStepControl ThisType;
171 typedef ImplicitRungeKuttaSolverParameters ParameterType;
173 ImplicitRungeKuttaTimeStepControl ( TimeProviderType &timeProvider,
const ParameterType ¶meters )
174 : timeProvider_( timeProvider ),
175 parameters_( parameters.clone() ),
176 cfl_( parameters_->cflStart() ),
177 cflMax_( parameters_->cflMax() ),
178 verbose_( parameters_->verbose() ),
179 initialized_( false )
182 explicit ImplicitRungeKuttaTimeStepControl ( TimeProviderType &timeProvider,
const Dune::Fem::ParameterReader ¶meter = 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 )
191 double time ()
const {
return timeProvider_.time(); }
192 double timeStepSize ()
const {
return timeProvider_.deltaT(); }
194 void initialTimeStepSize (
double helmholtzEstimate,
double sourceTermEstimate )
197 parameters().initTimeStepEstimate( sourceTermEstimate, helmholtzEstimate, estimate, cfl_ );
198 timeProvider_.provideTimeStepEstimate( estimate );
202 template<
class Monitor >
203 void reduceTimeStep (
double helmholtzEstimate,
double sourceTermEstimate,
const Monitor &monitor )
209 parameters().cflFactor( helmholtzEstimate, sourceTermEstimate, monitor.linearSolverIterations_,
false, factor );
216 if( (verbose_ >= ImplicitRungeKuttaSolverParameters::noConvergenceVerbosity) && (Dune::Fem::MPIManager::rank() == 0) )
218 Dune::derr <<
"Implicit Runge-Kutta step failed to converge." << std::endl;
220 <<
", iterations per time step: ILS = " << monitor.linearSolverIterations_
221 <<
", INLS = " << monitor.newtonIterations_
225 timeProvider_.provideTimeStepEstimate( factor * timeStepSize() );
226 timeProvider_.invalidateTimeStep();
229 template<
class Monitor >
230 void timeStepEstimate (
double helmholtzEstimate,
double sourceTermEstimate,
const Monitor &monitor )
238 parameters().cflFactor( helmholtzEstimate, sourceTermEstimate, monitor.linearSolverIterations_,
true, factor );
242 const double oldCfl = cfl_;
243 cfl_ = std::min( cflMax_, factor * cfl_ );
245 timeProvider_.provideTimeStepEstimate( std::min( sourceTermEstimate, cfl_ * helmholtzEstimate ) );
247 if( (cfl_ != oldCfl) && (verbose_ >= ImplicitRungeKuttaSolverParameters::cflVerbosity) && (Dune::Fem::MPIManager::rank() == 0) )
250 <<
", iterations per time step: ILS = " << monitor.linearSolverIterations_
251 <<
", INLS = " << monitor.newtonIterations_
256 bool computeError ()
const {
return false; }
259 const ParameterType ¶meters ()
const
261 assert( parameters_ );
265 TimeProviderType &timeProvider_;
266 std::shared_ptr< const ParameterType > parameters_;
267 double cfl_, cflMax_;
290 typedef ImplicitRungeKuttaTimeStepControl BaseType;
293 using BaseType :: initialized_;
294 using BaseType :: cfl_;
295 using BaseType :: parameters ;
298 typedef typename BaseType::ParameterType ParameterType;
301 : BaseType( timeProvider, parameters ),
305 if( parameters.parameter().getValue(
"fem.ode.pidcontrol",
bool(
false) ) )
307 tol_ = parameters.parameter().getValue(
"fem.ode.pidtolerance", tol_ );
308 errors_.resize( 3, tol_ );
313 : BaseType( timeProvider, parameter ),
317 if( parameter.getValue(
"fem.ode.pidcontrol",
bool(
false) ) )
319 tol_ = parameter.getValue(
"fem.ode.pidtolerance", tol_ );
320 errors_.resize( 3, tol_ );
324 bool computeError ()
const {
return ! errors_.empty() ; }
326 template<
class Monitor >
327 void timeStepEstimate (
double helmholtzEstimate,
double sourceTermEstimate,
const Monitor &monitor )
335 double dtEst = pidTimeStepControl( std::min( sourceTermEstimate, helmholtzEstimate ), monitor );
344 std::cout <<
"Set dt = " << dtEst << std::endl;
345 timeProvider_.provideTimeStepEstimate( dtEst );
347 if( (verbose_ >= ImplicitRungeKuttaSolverParameters::cflVerbosity) && (Dune::Fem::MPIManager::rank() == 0) )
350 <<
", iterations per time step: ILS = " << monitor.linearSolverIterations_
351 <<
", INLS = " << monitor.newtonIterations_
357 BaseType::timeStepEstimate( helmholtzEstimate, sourceTermEstimate, monitor );
361 template <
class Monitor >
362 double pidTimeStepControl(
const double dt,
const Monitor& monitor )
365 const double error = monitor.error_;
366 std::cout << error <<
" error " << std::endl;
367 if( std::abs( error ) < 1e-12 )
return 10. * dt;
370 for(
int i=0; i<2; ++i )
372 errors_[ i ] = errors_[i+1];
376 errors_[ 2 ] = error ;
381 const double newDt = dt * tol_ / error;
384 else if( error > 1e-12 )
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;
401 std::vector< double > errors_;
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