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...
|
| ERKStep ()=default |
| No memory allocation. More...
|
|
| ERKStep (ConvertsToButcherTableau< value_type > tableau, const ContainerType ©able) |
| 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...
|
|
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
◆ 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 >
the value type of the time variable (float or double)
◆ ERKStep() [1/2]
template<class ContainerType >
◆ ERKStep() [2/2]
template<class ContainerType >
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) |
◆ 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 >
All subsequent calls to step
method will enable the check for the first same as last property.
◆ ignore_fsal()
template<class ContainerType >
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 >
global order of the method given by the current Butcher Tableau
◆ step() [1/2]
template<class ContainerType >
template<class ExplicitRHS >
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 >
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: