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

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
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::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
tableauTableau, name or identifier that ConvertsToButcherTableau
copyablevector 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
Paramsdeduced by the compiler
Parameters
psparameters 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
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
rhsright hand side subroutine
t0start time
u0value 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)
dttimestep
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
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
rhsright hand side subroutine
t0start time
u0value 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)
dttimestep
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
deltaContains error estimate (u1 - tilde u1) on return (must have equal size as u0)

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