Embedded Runge Kutta explicit time-step with error estimate \( \begin{align} k_i = f\left( t^n + c_i \Delta t, u^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j\right) \\ u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s b_j k_j \\ \tilde u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s \tilde b_j k_j \\ \delta^{n+1} = u^{n+1} - \tilde u^{n+1} = \Delta t\sum_{j=1}^s (b_j - \tilde b_j) k_j \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
	ERKStep ()=default
	No memory allocation. More...

	ERKStep (ConvertsToButcherTableau< value_type > tableau, const ContainerType &copyable)
	Reserve internal workspace 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...

void	ignore_fsal ()
	All subsequent calls to `step` method will ignore the first same as last property (useful if you want to implement an operator splitting) More...

void	enable_fsal ()
	All subsequent calls to `step` method will enable the check for the first same as last property. More...

template<class ExplicitRHS >
void	step (ExplicitRHS &rhs, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt, ContainerType &delta)
	Advance one step with error estimate. More...

template<class ExplicitRHS >
void	step (ExplicitRHS &rhs, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt)
	Advance one step ignoring error estimate and embedded method. More...

unsigned	order () const
	global order of the method given by the current Butcher Tableau More...

unsigned	embedded_order () const
	global order of the embedding given by the current Butcher Tableau More...

unsigned	num_stages () const
	number of stages of the method given by the current Butcher Tableau More...

Detailed Description

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

Embedded Runge Kutta explicit time-step with error estimate \( \begin{align} k_i = f\left( t^n + c_i \Delta t, u^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j\right) \\ u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s b_j k_j \\ \tilde u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s \tilde b_j k_j \\ \delta^{n+1} = u^{n+1} - \tilde u^{n+1} = \Delta t\sum_{j=1}^s (b_j - \tilde b_j) k_j \end{align} \).

The method is defined by its (extended explicit) ButcherTableau, given by the coefficients a, b and c, and s is the number of stages. The embedding is given by the coefficients bt (tilde b).

You can provide your own coefficients or use one of the embedded methods in the following table:

We follow the naming convention of the ARKode library https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html (They also provide nice stability plots for their methods) as NAME-S-P-Q or NAME-S-Q, where

NAME is the author or name of the method
S is the number of stages in the method
P is the global order of the embedding
Q is the global order of the method

Name	Identifier	Description
Euler	dg::EXPLICIT_EULER_1_1	Explicit Euler
Midpoint-2-2	dg::MIDPOINT_2_2	Midpoint method
Kutta-3-3	dg::KUTTA_3_3	Kutta-3-3
Runge-Kutta-4-4	dg::CLASSIC_4_4	"The" Runge-Kutta method
SSPRK-2-2	dg::SSPRK_2_2	SSPRK (2,2) CFL_eff = 0.5 "Shu-Osher-Form"
SSPRK-3-2	dg::SSPRK_3_2	SSPRK (3,2) CFL_eff = 0.66 "Shu-Osher-Form"
SSPRK-3-3	dg::SSPRK_3_3	SSPRK (3,3) CFL_eff = 0.33 "Shu-Osher-Form"
SSPRK-5-3	dg::SSPRK_5_3	SSPRK (5,3) CFL_eff = 0.5 "Shu-Osher-Form"
SSPRK-5-4	dg::SSPRK_5_4	SSPRK (5,4) CFL_eff = 0.37 "Shu-Osher-Form"
Heun-Euler-2-1-2	dg::HEUN_EULER_2_1_2	Heun-Euler-2-1-2
Cavaglieri-3-1-2 (explicit)	dg::CAVAGLIERI_3_1_2	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems IMEXRKCB2 scheme
Fehlberg-3-2-3	dg::FEHLBERG_3_2_3	The original uses the embedding as the solution [Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987]
Fehlberg-4-2-3	dg::FEHLBERG_4_2_3	The original uses the embedding as the solution [Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987] (fsal)
Bogacki-Shampine-4-2-3	dg::BOGACKI_SHAMPINE_4_2_3	Bogacki-Shampine (fsal)
Cavaglieri-4-2-3 (explicit)	dg::CAVAGLIERI_4_2_3	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems The SSP scheme IMEXRKCB3c (explicit part)
ARK-4-2-3 (explicit)	dg::ARK324L2SA_ERK_4_2_3	ARK-4-2-3 (explicit)
Zonneveld-5-3-4	dg::ZONNEVELD_5_3_4	Zonneveld-5-3-4
ARK-6-3-4 (explicit)	dg::ARK436L2SA_ERK_6_3_4	ARK-6-3-4 (explicit)
Sayfy_Aburub-6-3-4	dg::SAYFY_ABURUB_6_3_4	Sayfy_Aburub_6_3_4
Cash_Karp-6-4-5	dg::CASH_KARP_6_4_5	Cash-Karp
Fehlberg-6-4-5	dg::FEHLBERG_6_4_5	Runge-Kutta-Fehlberg
Dormand-Prince-7-4-5	dg::DORMAND_PRINCE_7_4_5	Dormand-Prince method (fsal)
Tsitouras09-7-4-5	dg::TSITOURAS09_7_4_5	Tsitouras 5(4) method from 2009 (fsal), The default method in Julia
Tsitouras11-7-4-5	dg::TSITOURAS11_7_4_5	Tsitouras 5(4) method from 2011 (fsal) Further improves Tsitouras09 (Note that in the paper b-bt is printed instead of bt)
ARK-8-4-5 (explicit)	dg::ARK548L2SA_ERK_8_4_5	Kennedy and Carpenter (2019) Optimum ARK_2 method (explicit part)
Verner-9-5-6	dg::VERNER_9_5_6	Verner-9-5-6 (fsal)
Verner-10-6-7	dg::VERNER_10_6_7	Verner-10-6-7
Fehlberg-13-7-8	dg::FEHLBERG_13_7_8	Fehlberg-13-7-8
Dormand-Prince-13-7-8	dg::DORMAND_PRINCE_13_7_8	[Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987]
Feagin-17-8-10	dg::FEAGIN_17_8_10	Feagin (The RK10(8) method)

Note: The name of this class is in reminiscence of the ARKode library https://sundials.readthedocs.io/en/latest/arkode/index.html

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

the type of the vector class in use

◆ value_type

template<class ContainerType >

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

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

Constructor & Destructor Documentation

◆ ERKStep() [1/2]

template<class ContainerType >

dg::ERKStep< ContainerType >::ERKStep ( )

default

No memory allocation.

◆ ERKStep() [2/2]

template<class ContainerType >

dg::ERKStep< ContainerType >::ERKStep	(	ConvertsToButcherTableau< value_type >	tableau,
		const ContainerType &	copyable
	)

inline

Reserve internal workspace for the integration.

Parameters

tableau	Tableau, name or identifier that `ConvertsToButcherTableau`
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::ERKStep< 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::ERKStep< 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

◆ embedded_order()

template<class ContainerType >

unsigned dg::ERKStep< ContainerType >::embedded_order ( ) const

inline

global order of the embedding given by the current Butcher Tableau

◆ enable_fsal()

template<class ContainerType >

void dg::ERKStep< ContainerType >::enable_fsal ( )

inline

All subsequent calls to step method will enable the check for the first same as last property.

◆ ignore_fsal()

template<class ContainerType >

void dg::ERKStep< ContainerType >::ignore_fsal ( )

inline

All subsequent calls to step method will ignore the first same as last property (useful if you want to implement an operator splitting)

◆ num_stages()

template<class ContainerType >

unsigned dg::ERKStep< ContainerType >::num_stages ( ) const

inline

number of stages of the method given by the current Butcher Tableau

◆ order()

template<class ContainerType >

unsigned dg::ERKStep< ContainerType >::order ( ) const

inline

global order of the method given by the current Butcher Tableau

◆ step() [1/2]

template<class ContainerType >

template<class ExplicitRHS >

void dg::ERKStep< ContainerType >::step	(	ExplicitRHS &	rhs,
		value_type	t0,
		const ContainerType &	u0,
		value_type &	t1,
		ContainerType &	u1,
		value_type	dt
	)

inline

Advance one step ignoring error estimate and embedded method.

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.

Parameters

rhs	right hand side subroutine
t0	start time
u0	value at `t0`
t1	(write only) end time ( equals `t0+dt` on return, may alias `t0`)
u1	(write only) contains result on return (may alias u0)
dt	timestep

Note: on return rhs(t1, u1) will be the last call to rhs (this is useful if ExplicitRHS holds state, which is then updated to the current timestep); About the first same as last property (fsal): Some Butcher tableaus (e.g. Dormand-Prince or Bogacki-Shampine) have the property that the last value k_s of a timestep is the same as the first value k_0 of the next timestep. This means that we can save one call to the right hand side. This property is automatically activated if tableau.isFsal() returns true and t0 equals t1 of the last call to step. You can deactivate it by calling the ignore_fsal() method, which is useful for splitting methods but increases the number of rhs calls by 1.

Attention: On the rare occasion where you want to change the type of ExplicitRHS from one step to the next the fsal property is interpreted wrongly and will lead to wrong results. You will need to either re-construct the object or set the ignore_fsal property before the next step.

◆ step() [2/2]

template<class ContainerType >

template<class ExplicitRHS >

void dg::ERKStep< ContainerType >::step	(	ExplicitRHS &	rhs,
		value_type	t0,
		const ContainerType &	u0,
		value_type &	t1,
		ContainerType &	u1,
		value_type	dt,
		ContainerType &	delta
	)

inline

Advance one step with error estimate.

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.

Parameters

rhs	right hand side subroutine
t0	start time
u0	value at `t0`
t1	(write only) end time ( equals `t0+dt` on return, may alias `t0`)
u1	(write only) contains result on return (may alias u0)
dt	timestep

Note: on return rhs(t1, u1) will be the last call to rhs (this is useful if ExplicitRHS holds state, which is then updated to the current timestep); About the first same as last property (fsal): Some Butcher tableaus (e.g. Dormand-Prince or Bogacki-Shampine) have the property that the last value k_s of a timestep is the same as the first value k_0 of the next timestep. This means that we can save one call to the right hand side. This property is automatically activated if tableau.isFsal() returns true and t0 equals t1 of the last call to step. You can deactivate it by calling the ignore_fsal() method, which is useful for splitting methods but increases the number of rhs calls by 1.

Attention: On the rare occasion where you want to change the type of ExplicitRHS from one step to the next the fsal property is interpreted wrongly and will lead to wrong results. You will need to either re-construct the object or set the ignore_fsal property before the next step.

Parameters

delta Contains error estimate (u1 - tilde u1) on return (must have equal size as u0)

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

runge_kutta.h

Public Types

Public Member Functions

Detailed Description

Member Typedef Documentation

◆ container_type

◆ value_type

Constructor & Destructor Documentation

◆ ERKStep() [1/2]

◆ ERKStep() [2/2]

Member Function Documentation

◆ construct()

◆ copyable()

◆ embedded_order()

◆ enable_fsal()

◆ ignore_fsal()

◆ num_stages()

◆ order()

◆ step() [1/2]

◆ step() [2/2]