Semi-implicit multistep ODE integrator \( \begin{align} v^{n+1} = \sum_{q=0}^{s-1} a_q v^{n-q} + \Delta t\left[\left(\sum_{q=0}^{s-1}b_q \hat E(t^{n}-q\Delta t, v^{n-q}) + \sum_{q=1}^{s} c_q \hat I( t^n - q\Delta t, v^{n-q})\right) + c_0\hat I(t^{n}+\Delta t, v^{n+1})\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
	ImExMultistep ()=default
	No memory allocation. More...

	ImExMultistep (ConvertsToMultistepTableau< value_type > tableau, const ContainerType &copyable)
	Reserve memory for integration and construct Solver. 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 , class ImplicitRHS , class Solver >
void	init (const std::tuple< ExplicitRHS, ImplicitRHS, Solver > &ode, value_type t0, const ContainerType &u0, value_type dt)
	Initialize timestepper. Call before using the step function. More...

template<class ExplicitRHS , class ImplicitRHS , class Solver >
void	step (const std::tuple< ExplicitRHS, ImplicitRHS, Solver > &ode, value_type &t, ContainerType &u)
	Advance one timestep. More...

Detailed Description

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

Semi-implicit multistep ODE integrator \( \begin{align} v^{n+1} = \sum_{q=0}^{s-1} a_q v^{n-q} + \Delta t\left[\left(\sum_{q=0}^{s-1}b_q \hat E(t^{n}-q\Delta t, v^{n-q}) + \sum_{q=1}^{s} c_q \hat I( t^n - q\Delta t, v^{n-q})\right) + c_0\hat I(t^{n}+\Delta t, v^{n+1})\right] \end{align} \).

which discretizes

\[ \frac{\partial v}{\partial t} = \hat E(t,v) + \hat I(t,v) \]

where \( \hat E \) contains the explicit and \( \hat I \) the implicit part of the equations. As an example, the coefficients for the 3-step, 3rd order "Karniadakis" scheme are

\[ 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} \\ c_0 = \frac{6}{11}\quad c_1 = c_2 = c_3 = 0 \\ \]

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
ImEx-Euler-1-1	dg::IMEX_EULER_1_1	Explicit Euler combined with Implicit Euler
Euler	dg::IMEX_EULER_1_1	For convenience
ImEx-Koto-2-2	dg::IMEX_KOTO_2_2	Koto T. Front. Math. China 2009, 4(1): 113-129 A stabilized 2nd order scheme with a large region of stability
ImEx-Adams-X-X	dg::IMEX_ADAMS_X_X	Hundsdorfer and Ruuth, Journal of Computational Physics 225 (2007) Note Possible values for X: 2 (C=0.44), 3 (C=0.16)
ImEx-BDF-X-X	dg::IMEX_BDF_X_X	The family of schems described in Hundsdorfer and Ruuth, Journal of Computational Physics 225 (2007) The implicit part is a normal BDF scheme https://en.wikipedia.org/wiki/Backward_differentiation_formula while the explicit part equals the Minimal Projecting method by Alfeld, P., Math. Comput. 33.148 1195-1212 (1979) or extrapolated BDF in Hundsdorfer, W., Ruuth, S. J., & Spiteri, R. J. (2003). Monotonicity-preserving linear multistep methods. SIAM Journal on Numerical Analysis, 41(2), 605-623 Note Possible values for X: 1 (C=1.00), 2 (C=0.63), 3 (C=0.39), 4 (C=0.22), 5 (C=0.09), 6 Note that X=3 is identical to the "Karniadakis" scheme
Karniadakis	dg::IMEX_BDF_3_3	The ImEx-BDF-3-3 scheme is identical to the widely used "Karniadakis" scheme Karniadakis, et al. J. Comput. Phys. 97 (1991)
ImEx-TVB-X-X	dg::IMEX_TVB_X_X	The family of schems described in < Hundsdorfer and Ruuth, Journal of Computational Physics 225 (2007) The explicit part is a TVB scheme while the implicit part is optimized to maximize damping of high wavelength Note Possible values for X: 3 (C=0.54), 4 (C=0.46), 5 (C=0.38)

Note: 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

The necessary Inversion in the implicit part is provided by the Implicit class.

The following code example demonstrates how to implement the method of manufactured solutions on a 2d partial differential equation with the dg library:

struct Explicit
{
    Explicit( double nu): m_nu( nu) {}
    void operator()( double t, const std::array<double,2>& T,
            std::array<double,2>& Tp)
    {
        Tp[0] = m_nu*cos(t) - sin(t);
        Tp[1] = m_nu*sin(t) + cos(t);
    }
    private:
    double m_nu;
};
 
//the implicit part contains  Tp = -nu T
struct Implicit
{
    Implicit( double nu): m_nu(nu) { }
 
    void operator()( double t, const std::array<double,2>& T,
            std::array<double,2>& Tp)
    {
        Tp[0] = -m_nu*T[0];
        Tp[1] = -m_nu*T[1];
    }
    void operator()( double alpha, double t, std::array<double,2>& y, const
            std::array<double,2>& rhs)
    {
        // solve y - alpha I(t,y)  = rhs
        y[0] = rhs[0]/(1+alpha*m_nu);
        y[1] = rhs[1]/(1+alpha*m_nu);
    }
  private:
    double m_nu;
};

In the main function:

        //construct time stepper
        dg::ImExMultistep< std::array<double,2> > imex( name, y0);
        time = 0., y0 = init; //y0 and init are of type std::array<double,2> and contain the initial condition
        //initialize the timestepper (ex and im are objects of type Explicit and Implicit defined above)
        imex.init( std::tie(ex, im, im), time, y0, dt);
        //main time loop (NT = 20)
        for( unsigned k=0; k<NT; k++)
            imex.step( std::tie(ex, im, im), time, y0); //inplace step

Note: In our experience the implicit treatment of diffusive or hyperdiffusive terms may significantly reduce the required number of time steps. This outweighs the increased computational cost of the additional matrix inversions. However, each PDE is different and general statements like this one should be treated with care.; 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

ContainerType

Any 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::ImExMultistep< ContainerType >::container_type = ContainerType

the type of the vector class in use

◆ value_type

template<class ContainerType >

using dg::ImExMultistep< ContainerType >::value_type = get_value_type<ContainerType>

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

Constructor & Destructor Documentation

◆ ImExMultistep() [1/2]

template<class ContainerType >

dg::ImExMultistep< ContainerType >::ImExMultistep ( )

default

No memory allocation.

◆ ImExMultistep() [2/2]

template<class ContainerType >

dg::ImExMultistep< ContainerType >::ImExMultistep	(	ConvertsToMultistepTableau< value_type >	tableau,
		const ContainerType &	copyable
	)

inline

Reserve memory for integration and construct Solver.

Parameters

tableau	Tableau, name or identifier that `ConvertsToMultistepTableau`
copyable	vector 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 the same size)

Member Function Documentation

◆ construct()

template<class ContainerType >

template<class ... Params>

void dg::ImExMultistep< ContainerType >::construct ( Params &&... ps )

inline

Perfect forward parameters to one of the constructors.

Template Parameters

Params deduced by the compiler

Parameters

ps	parameters forwarded to constructors

◆ copyable()

template<class ContainerType >

const ContainerType & dg::ImExMultistep< 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 , class ImplicitRHS , class Solver >

void dg::ImExMultistep< ContainerType >::init	(	const std::tuple< ExplicitRHS, ImplicitRHS, Solver > &	ode,
		value_type	t0,
		const ContainerType &	u0,
		value_type	dt
	)

Initialize timestepper. Call before using the step function.

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

Template Parameters

ExplicitRHS	The 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`.
ImplicitRHS	The implicit (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' = I(t, y) translates to I(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`.
Solver	A functor type for the implicit right hand side. Must solve the equation \( y - \alpha I(y,t) = y^*\) for `y` for given `alpha`, `t` and `ys`. Alpha is always positive and non-zero. Signature `void operator()( value_type alpha, value_type t, ContainerType& y, const ContainerType& ys);` The functor can throw. Any Exception should derive from `std::exception`.

Parameters

ode	the <explicit rhs, implicit rhs, solver for the rhs> functor. Typically `std::tie(explicit_rhs, implicit_rhs, solver)`
t0	The intital time corresponding to u0
u0	The initial value of the integration
dt	The timestep saved for later use

Note: the implementation is such that on return the last call is the explicit part explicit_part at (t0,u0). This is useful if explicit_part holds state, which is then updated to that timestep and/or if implicit_rhs changes the state of explicit_rhs

◆ step()

template<class ContainerType >

template<class ExplicitRHS , class ImplicitRHS , class Solver >

void dg::ImExMultistep< ContainerType >::step	(	const std::tuple< ExplicitRHS, ImplicitRHS, Solver > &	ode,
		value_type &	t,
		ContainerType &	u
	)

Advance one timestep.

Template Parameters

ExplicitRHS	The 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`.
ImplicitRHS	The implicit (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' = I(t, y) translates to I(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`.
Solver	A functor type for the implicit right hand side. Must solve the equation \( y - \alpha I(y,t) = y^*\) for `y` for given `alpha`, `t` and `ys`. Alpha is always positive and non-zero. Signature `void operator()( value_type alpha, value_type t, ContainerType& y, const ContainerType& ys);` The functor can throw. Any Exception should derive from `std::exception`.

Parameters

ode	the <explicit rhs, implicit rhs, solver for the rhs> functor. Typically `std::tie(explicit_rhs, implicit_rhs, solver)`
t	(write-only), contains timestep corresponding to `u` on return
u	(write-only), contains next step of time-integration on return

Note: After a solve, we call explicit_rhs at the new (t,u) which is the last call before return. This is useful if explicit_rhs holds state, which is then updated to the new timestep and/or if solver changes the state of explicit_rhs

Attention: The first few steps after the call to the init function are performed with a semi-implicit Runge-Kutta method to initialize the multistepper

Note: The only time the implicit_rhs functor is called is during the initialization phase (the first few steps after the call to the init function)

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

multistep.h

Public Types

Public Member Functions

Detailed Description

Member Typedef Documentation

◆ container_type

◆ value_type

Constructor & Destructor Documentation

◆ ImExMultistep() [1/2]

◆ ImExMultistep() [2/2]

Member Function Documentation

◆ construct()

◆ copyable()

◆ init()

◆ step()