|
| | Adaptive ()=default |
| |
| template<class ... StepperParams> |
| | Adaptive (StepperParams &&...ps) |
| | Allocate workspace and construct stepper. More...
|
| |
| template<class ... Params> |
| void | construct (Params &&...ps) |
| | Perfect forward parameters to one of the constructors. More...
|
| |
| stepper_type & | stepper () |
| | Allow write access to internal stepper. More...
|
| |
| const stepper_type & | stepper () const |
| | Read access to internal stepper. More...
|
| |
| template<class ODE , class ControlFunction = value_type (std::array<value_type,3>, std::array< value_type,3>, unsigned , unsigned), class ErrorNorm = value_type( const container_type&)> |
| void | step (ODE &&ode, value_type t0, const container_type &u0, value_type &t1, container_type &u1, value_type &dt, ControlFunction control, ErrorNorm norm, value_type rtol, value_type atol, value_type reject_limit=2) |
| | Semi-implicit adaptive step. More...
|
| |
| bool | failed () const |
| | Return true if the last stepsize in step was rejected. More...
|
| |
| const unsigned & | nfailed () const |
| | Get total number of failed steps. More...
|
| |
| unsigned & | nfailed () |
| | Set total number of failed steps. More...
|
| |
| const unsigned & | nsteps () const |
| | Get total number of step calls. More...
|
| |
| unsigned & | nsteps () |
| | Re-set total number of step calls. More...
|
| |
| const value_type & | get_error () const |
| | Get the latest error norm relative to solution vector. More...
|
| |
template<class Stepper>
struct dg::Adaptive< Stepper >
Driver class for adaptive timestep ODE integration.
In order to build an adaptive ODE integrator you basically need three ingredients: a Stepper, a ControlFunction and an ErrorNorm. The Stepper does the actual computation and advances the solution one step further with a given timestep dt. Furthermore, it has to come up with an estimate of the error of the solution \(\delta_{n+1}\) and indicate the order of that error. With the ErrorNorm the error estimate can be converted to a scalar that can be compared to given relative and absolute error tolerances rtol and atol.
\[ \epsilon_{n+1} = || \frac{\delta_{n+1}}{(\epsilon_{rtol} |u_{n}| + \epsilon_{atol})\sqrt{N}}|| \]
where N is the array size, n is the solution at the previous timestep and the fraction is to be understood as a pointwise division of the vector elements. The ControlFunction will try to keep \( \epsilon_{n+1}\) close to 1 and comes up with an adapted suggestion for the timestep in the next step.
This form of error control entails that on average every point in the solution vector on its own fulfills \( |\delta_{n+1,i}| \approx \epsilon_{rtol}|u_{n,i} + \epsilon_{atol}\), (which places emphasize on atol in regions where the solution is close to zero).
However, if \(\epsilon_{n+1} > r\) where r=2 by default is the user-adaptable reject-limit, the step is rejected and the step will be recomputed and the controller restarted. For more information on these concepts we recommend the mathematical primer of the ARKode library.
For an example on how to use this class in a practical example consider the following code snippet:
time = 0., y0 = init;
double time = 0;
double dt = 1e-6;
int counter=0;
while( time < T )
{
if( time + dt > T)
dt = T-time;
adapt.step( std::tie(ex, im, im), time, y0, time, y0, dt,
counter ++;
}
static auto imex_control
h_{n+1} = |ex_control| < |im_control| ? ex_control : im_control
Definition: adaptive.h:109
static auto l2norm
Compute using dg::blas1::dot.
Definition: adaptive.h:22
Driver class for adaptive timestep ODE integration.
Definition: adaptive.h:233
- Template Parameters
-
| Stepper | A timestepper class that computes the actual timestep and an error estimate, for example an embedded Runge Kutta method dg::ERKStep or the additive method dg::ARKStep. But really, you can also choose to use your own timestepper class. The requirement is that there is a step member function that is called as stepper.step( ode, t0, u0, t1, u1, dt, delta) Here, t0, t1 and dt are of type Stepper::value_type, u0,u1 and delta are vector types of type Stepper::container_type& and ode, ex and im are functors implementing the equations that are forwarded from the caller. The parameters t1, u1 and delta are output parameters and must be updated by the stepper. The Stepper must have the order() and embedded_order() member functions that return the (global) order of the method and its error estimate. The const ContainerType& copyable()const; member must return a container of the size that is later used in step (it does not matter what values copyable contains, but its size is important; the step method can only be called with vectors of this size) |
- Note
- On step rejection, choosing timesteps and introducing restrictions on the controller: here is a quote from professor G. Söderlind (the master of control functions) from a private e-mail conversation:
-
"The issue is that most controllers are best left to do their work
without interference. Each time one interferes with the control loop, one
creates a transient that has to be dealt with by the controller.
@note It is indeed necessary to reject steps every once in a while, but I
usually try to avoid it to the greatest possible extent. In my opinion, the
fear of having a too large error on some steps is vastly exaggerated. You see,
in some steps the error is too large, and in others it is too small, and all
controllers I constructed are designed to be “expectation value correct” in the
sense that if the errors are random, the too large and too small errors
basically cancel in the long run.
@note Even so, there are times when the error is way out of proportion. But I
usually accept an error that is up to, say 2*TOL, which typically won’t cause
any problems. Of course, if one hits a sharp change in the solution, the error
may be much larger, and the step recomputed. But then one must “reset" the controller, i.e., the controller keeps back information, and it is better to restart the controller to avoid back information to affect the next step."
- Attention
- Should you use this class instead of a fixed stepsize Multistep say? The thing to consider especially when solving partial differential equations, is that the right hand side might be very costly to evaluate. An adaptive stepper (especially one-step methods) usually calls this right hand side more often than a multistep (only one call per step). The additional computation of the error norm in the Adaptive step might also be important since the norm requires global communication in a parallel program (bad for scaling to many nodes). So does the Adaptive timestepper make up for its increased cost throught the adaption of the stepsize? In some cases the answer might be a sound Yes. Especially when there are velocity peaks in the solution the multistep timestep might be restricted too much. In other cases when the timestep does not need to be adapted much, a multistep method can be faster. In any case the big advantage of Adaptive is that it usually just works (even though it is not fool-proof) and you do not have to spend time finding a suitable timestep like in the multistep method.
template<class Stepper >
template<class ODE , class ControlFunction = value_type (std::array<value_type,3>, std::array< value_type,3>, unsigned , unsigned), class ErrorNorm = value_type( const container_type&)>
| void dg::Adaptive< Stepper >::step |
( |
ODE && |
ode, |
|
|
value_type |
t0, |
|
|
const container_type & |
u0, |
|
|
value_type & |
t1, |
|
|
container_type & |
u1, |
|
|
value_type & |
dt, |
|
|
ControlFunction |
control, |
|
|
ErrorNorm |
norm, |
|
|
value_type |
rtol, |
|
|
value_type |
atol, |
|
|
value_type |
reject_limit = 2 |
|
) |
| |
|
inline |
Semi-implicit adaptive step.
- Parameters
-
| t0 | initial time |
| u0 | initial value at t0 |
| t1 | (write only) end time ( equals t0+dt on output if the step was accepted, otherwise equals t0, may alias t0) |
| u1 | (write only) contains the updated result on output if the step was accepted, otherwise a copy of u0 (may alias u0) |
| dt | on input: timestep On output: stepsize proposed by the controller that can be used to continue the integration in the next step. |
- Template Parameters
-
| ODE | The ExplicitRHS or tuple type that corresponds to what is inserted into the step member of the Stepper |
- Parameters
-
- Template Parameters
-
| ControlFunction | function or Functor called as dt' = control( {dt0, dt1, dt2}, {eps0, eps1, eps2}, embedded_order, order), where all parameters are of type value_type except the last two, which are unsigned. |
- Parameters
-
| control | The control function. For explicit and imex methods, dg::pid_control is a good choice with dg::ex_control or dg::imex_control as an alternative if too many steps fail. For implicit methods use the dg::im_control. The task of the control function is to compute a new timestep size based on the old timestep size, the order of the method and the past error(s). The behaviour of the controller is also changed by the set_reject_limit function |
- Template Parameters
-
| ErrorNorm | function or Functor of type value_type( const
ContainerType&) |
- Parameters
-
| norm | The error norm. Usually dg::l2norm is a good choice, but for very small vector sizes the time for the binary reproducible dot product might become a performance bottleneck. Then dg::fast_l2norm is a better choice. |
| rtol | the desired relative accuracy. Usually 1e-5 is a good choice. |
| atol | the desired absolute accuracy. Usually 1e-7 is a good choice. |
| reject_limit | the default value is 2. Sometimes even 2 is not enough, for example when the stepsize fluctuates very much and the stepper fails often then it may help to increase the reject_limit further (to 10 say). |
- Note
- The error tolerance is computed such that on average every point in the solution vector fulfills \( |\delta_i| \approx r|u_i| + a\)
-
Try not to mess with dt. The controller is best left alone and it does a very good job choosing timesteps. But how do I output my solution at certain (equidistant) timesteps? First, think about if you really, really need that. Why is it so bad to have output at non-equidistant timesteps? If you are still firm, then consider using an interpolation scheme (cf.
dg::Extrapolation). Let choosing the timestep yourself be the very last option if the others are not viable
-
From Kennedy and Carpenter, Appl. num. Math., (2003): "Step-size control is a means by which accuracy, iteration, and to a lesser extent stability are controlled.
The choice of (t) may be chosen from many criteria, among those are the (t) from the accuracy
based step controller, the (t)_inviscid and (t)_viscous associated with the inviscid and viscous stability
limits of the ERK, and the (t)_iter associated with iteration convergence."
Our control method is an error based control method. However, due to the CFL (or other stability) condition there might be a sharp barrier in the range of possible stepsizes, i.e the stability limited timestep where the error sharply increses leading to a large number of rejected steps. The pid-controller usually does a good job keeping the timestep "just right" but in future work we may want to consider stability based control methods as well. Furthermore, for implicit or semi-implicit methods the chosen timestep influences how fast the iterative solver converges. If the timestep is too large the solver might take too long to converge and a smaller timestep might be preferable. Even for explicit methods where in each timestep an elliptic equation has to be solved the timestep may influence the convergence rate. Currently, there is no communication implemented between these solvers and the controller, but we may consider it in future releases.
- Attention
- When you use the ARKStep in combination with the Adaptive time step algorithm pay attention to solve the implicit part with sufficient accuracy. Else, the error propagates into the time controller, which will then choose the timestep as if the implicit part was explicit i.e. far too small. This might have to do with stiffness-leakage [Kennedy and Carpenter, Appl. num. Math., (2003)]: "An essential requirement for the viability of stiff/nonstiff IMEX schemes is that the stiffness remains
truely separable. If this were not the case then stiffness would leak out of the stiff terms and stiffen the
nonstiff terms. It would manifest itself as a loss in stability or a forced reduction in stepsize of the nonstiff
terms. A more expensive fully implicit approach might then be required, and hence, methods that leak
substantial stiffness might best be avoided".
- Note
- The exact values of
rtol and atol might not be too important. However, don't make rtol too large, 1e-1 say, since then the controller might get too close to the CFL barrier. The timestepper is still able to crash, mind, even though the chances of that happening are somewhat lower than in a fixed stepsize method.
1.9.3