Discontinuous Galerkin Library
#include "dg/algorithm.h"
dg::ExplicitMultistep< ContainerType > Struct Template Reference

General explicit linear multistep ODE integrator \( \begin{align} v^{n+1} = \sum_{j=0}^{s-1} a_j v^{n-j} + \Delta t\left(\sum_{j=0}^{s-1}b_j \hat f\left(t^{n}-j\Delta t, v^{n-j}\right)\right) \end{align} \). More...

Public Types

using value_type = get_value_type< ContainerType >
 the value type of the time variable (float or double) More...
 
using container_type = ContainerType
 

Public Member Functions

 ExplicitMultistep ()=default
 No memory allocation. More...
 
 ExplicitMultistep (ConvertsToMultistepTableau< value_type > tableau, const ContainerType &copyable)
 Reserve memory for the integration. More...
 
template<class ... Params>
void construct (Params &&...ps)
 Perfect forward parameters to one of the constructors. More...
 
const ContainerType & copyable () const
 Return an object of same size as the object used for construction. More...
 
template<class ExplicitRHS >
void init (ExplicitRHS &rhs, value_type t0, const ContainerType &u0, value_type dt)
 Initialize timestepper. Call before using the step function. More...
 
template<class ExplicitRHS >
void step (ExplicitRHS &rhs, value_type &t, ContainerType &u)
 Advance one timestep. More...
 

Detailed Description

template<class ContainerType>
struct dg::ExplicitMultistep< ContainerType >

General explicit linear multistep ODE integrator \( \begin{align} v^{n+1} = \sum_{j=0}^{s-1} a_j v^{n-j} + \Delta t\left(\sum_{j=0}^{s-1}b_j \hat f\left(t^{n}-j\Delta t, v^{n-j}\right)\right) \end{align} \).

which discretizes

\[ \frac{\partial v}{\partial t} = \hat f(t,v) \]

where \( f \) contains the equations. The coefficients for an order 3 "eBDF" scheme are given as an example:

\[ a_0 = \frac{18}{11}\ a_1 = -\frac{9}{11}\ a_2 = \frac{2}{11} \\ b_0 = \frac{18}{11}\ b_1 = -\frac{18}{11}\ b_2 = \frac{6}{11} \]

You can use your own coefficients defined as a dg::MultistepTableau or use one of the predefined coefficients in

We follow the naming convention as NAME-S-Q

  • NAME is the author or name of the method
  • S is the number of steps in the method
  • Q is the global order of the method
Name Identifier Description
AB-X-X dg::AB_X_X The family of schemes described in Linear multistep methods as Adams-Bashforth

\[ u^{n+1} = u^n + \Delta t\sum_{j=0}^{s-1} b_j f\left(t^n - j \Delta t, u^{n-j}\right) \]

Note
Possible stages are X: 1, 2,..., 5, the order of the method is the same as its stages
The Adams-Bashforth schemes implemented here need less storage but may have a smaller region of absolute stability than for example an extrapolated BDF method of the same order.
eBDF-X-X dg::eBDF_X_X The family of schemes described in Hundsdorfer, W., Ruuth, S. J., & Spiteri, R. J. (2003). Monotonicity-preserving linear multistep methods. SIAM Journal on Numerical Analysis, 41(2), 605-623 as extrapolated BDF where it is found to be TVB (total variation bound). The schemes also appear as Minimal Projecting scheme described in Alfeld, P., Math. Comput. 33.148 1195-1212 (1979)
Note
Possible stages are X: 1 (C=1), 2 (C=0.63), 3 (C=0.39), 4 (C=0.22), 5 (C=0.09), 6 with the order the same as the number of stages
TVB-X-X dg::TVB_X_X The family of schemes described in S.J. Ruuth and W. Hundsdorfer, High-order linear multistep methods with general monotonicity and boundedness properties, Journal of Computational Physics, Volume 209, Issue 1, 2005 as Total variation Bound. These schemes have larger allowable step sizes than the eBDF family,
Note
Possible values for X are 1 (C=1), 2 (C=0.5), 3 (C=0.54), 4 (C=0.46), 5 (C=0.38) 6 (C=0.33). We highlight that TVB-3-3 has 38% larger allowable stepsize than eBDF-3-3 and TVB-4-4 has 109% larger stepsize than eBDF-4-4 (to ensure the TVB property, not stability).
SSP-X-Y dg::SSP_X_Y The family of schemes described in Gottlieb, S. On high order strong stability preserving runge-kutta and multi step time discretizations. J Sci Comput 25, 105–128 (2005) as Strong Stability preserving. We implement the lowest order schemes for each stage and disregard the remaining schemes in the paper since their CFL conditions are worse than the TVB scheme of the same order.
Note
Possible values for X-Y : 1-1 (C=1), 2-2 (C=0.5), 3-2 (C=0.5), 4-2 (C=0.66), 5-3 (C=0.5), 6-3 (C=0.567).
These schemes are noteworthy because the coefficients b_i are all positive except for the 2-2 method and the "4-2" and "6-3" methods allow slightly larger allowable stepsize but increased storage requirements than TVB of same order (2 and 3).
Note
Total variation bound (TVB) means \( || v^n|| \leq M ||v^0||\) where the norm signifies the total variation semi-norm. Total variation diminishing (TVD) means M=1, and strong stability preserving (SSP) is the same as TVD, TVB schemes converge to the correct entropy solutions of hyperbolic conservation laws
the CFL coefficient C is given relative to the forward Euler method: \( \Delta t < C \Delta t_{FE}\).
Attention
The coefficient C is the one that ensures the TVD property of the scheme and is not directly related to the stability region of the scheme
Note
Uses only blas1::axpby routines to integrate one step.
The difference between a multistep and a single step method like RungeKutta is that the multistep only takes one right-hand-side evaluation per step. This is advantageous if the right hand side is expensive to evaluate. Even though Runge Kutta methods can have a larger absolute timestep, if the effective timestep per rhs evaluation is compared, multistep methods generally win.
a disadvantage of multistep is that timestep adaption is not easily done.
Template Parameters
ContainerTypeAny class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others
  • dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp)
  • std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors)
  • double (scalar) and other primitive types ...
If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them
See also
See The dg dispatch system for a detailed explanation of our type dispatch system

Member Typedef Documentation

◆ container_type

template<class ContainerType >
using dg::ExplicitMultistep< ContainerType >::container_type = ContainerType

the type of the vector class in use

◆ value_type

template<class ContainerType >
using dg::ExplicitMultistep< ContainerType >::value_type = get_value_type<ContainerType>

the value type of the time variable (float or double)

Constructor & Destructor Documentation

◆ ExplicitMultistep() [1/2]

template<class ContainerType >
dg::ExplicitMultistep< ContainerType >::ExplicitMultistep ( )
default

No memory allocation.

◆ ExplicitMultistep() [2/2]

template<class ContainerType >
dg::ExplicitMultistep< ContainerType >::ExplicitMultistep ( ConvertsToMultistepTableau< value_type tableau,
const ContainerType &  copyable 
)
inline

Reserve memory for the integration.

Set the coefficients \( a_i,\ b_i\)

Parameters
tableauTableau, name or identifier that ConvertsToMultistepTableau
copyableContainerType of the size that is used in step
Note
it does not matter what values copyable contains, but its size is important

Member Function Documentation

◆ construct()

template<class ContainerType >
template<class ... Params>
void dg::ExplicitMultistep< ContainerType >::construct ( Params &&...  ps)
inline

Perfect forward parameters to one of the constructors.

Template Parameters
Paramsdeduced by the compiler
Parameters
psparameters forwarded to constructors

◆ copyable()

template<class ContainerType >
const ContainerType & dg::ExplicitMultistep< ContainerType >::copyable ( ) const
inline

Return an object of same size as the object used for construction.

Returns
A copyable object; what it contains is undefined, its size is important

◆ init()

template<class ContainerType >
template<class ExplicitRHS >
void dg::ExplicitMultistep< ContainerType >::init ( ExplicitRHS &  rhs,
value_type  t0,
const ContainerType &  u0,
value_type  dt 
)
inline

Initialize timestepper. Call before using the step function.

This routine has to be called before the first timestep is made.

Template Parameters
ExplicitRHSThe explicit (part of the) right hand side is a functor type with no return value (subroutine) of signature void operator()(value_type, const ContainerType&, ContainerType&) The first argument is the time, the second is the input vector, which the functor may not override, and the third is the output, i.e. y' = E(t, y) translates to E(t, y, y'). The two ContainerType arguments never alias each other in calls to the functor. The functor can throw to indicate failure. Exceptions should derive from std::exception.
Parameters
rhsThe rhs functor
t0The intital time corresponding to u0
u0The initial value of the integration
dtThe timestep saved for later use
Note
the implementation is such that on return the last call to the explicit part ex is at (t0,u0). This might be interesting if the call to ex changes its state.

◆ step()

template<class ContainerType >
template<class ExplicitRHS >
void dg::ExplicitMultistep< ContainerType >::step ( ExplicitRHS &  rhs,
value_type t,
ContainerType &  u 
)
inline

Advance one timestep.

Template Parameters
ExplicitRHSThe explicit (part of the) right hand side is a functor type with no return value (subroutine) of signature void operator()(value_type, const ContainerType&, ContainerType&) The first argument is the time, the second is the input vector, which the functor may not override, and the third is the output, i.e. y' = E(t, y) translates to E(t, y, y'). The two ContainerType arguments never alias each other in calls to the functor. The functor can throw to indicate failure. Exceptions should derive from std::exception.
Parameters
rhsThe rhs functor
t(write-only), contains timestep corresponding to u on return
u(write-only), contains next step of time-integration on return
Note
the implementation is such that on return the last call to the explicit part ex is at the new (t,u). This might be interesting if the call to ex changes its state.
Attention
The first few steps after the call to the init function are performed with a Runge-Kutta method (of the same order) to initialize the multistepper

The documentation for this struct was generated from the following file: