Utilities for ODE integrators

Collaboration diagram for Utilities for ODE integrators:

Classes
struct	dg::EntireDomain
	The domain that contains all points. More...
struct	dg::AdaptiveTimeloop< ContainerType >
	Integrate using a while loop. More...
struct	dg::MultistepTimeloop< ContainerType >
	Integrate using a for loop and a fixed non-changeable time-step. More...
struct	dg::MultistepTableau< real_type >
	Manage coefficients of Multistep methods. More...
struct	dg::ConvertsToMultistepTableau< real_type >
	Convert identifiers to their corresponding dg::MultistepTableau. More...
struct	dg::aTimeloop< ContainerType >
	Abstract timeloop independent of stepper and ODE. More...
struct	dg::IdentityFilter
	A filter that does nothing. More...
struct	dg::SinglestepTimeloop< ContainerType >
	Integrate using a for loop and a fixed time-step. More...
struct	dg::ConvertsToButcherTableau< real_type >
	Convert identifiers to their corresponding dg::ButcherTableau. More...
struct	dg::ConvertsToShuOsherTableau< real_type >
	Convert identifiers to their corresponding dg::ShuOsherTableau. More...

Enumerations
enum	dg::multistep_identifier {   dg::IMEX_EULER_1_1 , dg::IMEX_ADAMS_2_2 , dg::IMEX_ADAMS_3_3 , dg::IMEX_KOTO_2_2 ,   dg::IMEX_BDF_2_2 , dg::IMEX_BDF_3_3 , dg::IMEX_BDF_4_4 , dg::IMEX_BDF_5_5 ,   dg::IMEX_BDF_6_6 , dg::IMEX_TVB_3_3 , dg::IMEX_TVB_4_4 , dg::IMEX_TVB_5_5 ,   dg::AB_1_1 , dg::AB_2_2 , dg::AB_3_3 , dg::AB_4_4 ,   dg::AB_5_5 , dg::eBDF_1_1 , dg::eBDF_2_2 , dg::eBDF_3_3 ,   dg::eBDF_4_4 , dg::eBDF_5_5 , dg::eBDF_6_6 , dg::TVB_1_1 ,   dg::TVB_2_2 , dg::TVB_3_3 , dg::TVB_4_4 , dg::TVB_5_5 ,   dg::TVB_6_6 , dg::SSP_1_1 , dg::SSP_2_2 , dg::SSP_3_2 ,   dg::SSP_4_2 , dg::SSP_5_3 , dg::SSP_6_3 , dg::BDF_1_1 ,   dg::BDF_2_2 , dg::BDF_3_3 , dg::BDF_4_4 , dg::BDF_5_5 ,   dg::BDF_6_6 }
	Identifiers for Multistep Tableaus. More...
enum class	dg::to { dg::to::exact , dg::to::at_least }
	Switch for the Timeloop integrate function. More...
enum	dg::tableau_identifier {   dg::EXPLICIT_EULER_1_1 , dg::MIDPOINT_2_2 , dg::KUTTA_3_3 , dg::CLASSIC_4_4 ,   dg::HEUN_EULER_2_1_2 , dg::CAVAGLIERI_3_1_2 , dg::FEHLBERG_3_2_3 , dg::FEHLBERG_4_2_3 ,   dg::BOGACKI_SHAMPINE_4_2_3 , dg::CAVAGLIERI_4_2_3 , dg::ARK324L2SA_ERK_4_2_3 , dg::ZONNEVELD_5_3_4 ,   dg::ARK436L2SA_ERK_6_3_4 , dg::SAYFY_ABURUB_6_3_4 , dg::CASH_KARP_6_4_5 , dg::FEHLBERG_6_4_5 ,   dg::DORMAND_PRINCE_7_4_5 , dg::TSITOURAS09_7_4_5 , dg::TSITOURAS11_7_4_5 , dg::ARK548L2SA_ERK_8_4_5 ,   dg::VERNER_9_5_6 , dg::VERNER_10_6_7 , dg::FEHLBERG_13_7_8 , dg::DORMAND_PRINCE_13_7_8 ,   dg::FEAGIN_17_8_10 , dg::IMPLICIT_EULER_1_1 , dg::IMPLICIT_MIDPOINT_1_2 , dg::TRAPEZOIDAL_2_2 ,   dg::SDIRK_2_1_2 , dg::CAVAGLIERI_IMPLICIT_3_1_2 , dg::BILLINGTON_3_3_2 , dg::TRBDF2_3_3_2 ,   dg::SANCHEZ_3_3 , dg::KVAERNO_4_2_3 , dg::SDIRK_4_2_3 , dg::CAVAGLIERI_IMPLICIT_4_2_3 ,   dg::ARK324L2SA_DIRK_4_2_3 , dg::SANCHEZ_3_4 , dg::CASH_5_2_4 , dg::CASH_5_3_4 ,   dg::SDIRK_5_3_4 , dg::KVAERNO_5_3_4 , dg::ARK436L2SA_DIRK_6_3_4 , dg::KVAERNO_7_4_5 ,   dg::SANCHEZ_6_5 , dg::ARK548L2SA_DIRK_8_4_5 , dg::SANCHEZ_7_6 , dg::SSPRK_2_2 ,   dg::SSPRK_3_2 , dg::SSPRK_3_3 , dg::SSPRK_5_3 , dg::SSPRK_5_4 }
	Identifiers for Butcher Tableaus. More...

Functions
static std::string	dg::to2str (enum to mode)
	Convert integration mode to string. More...

Variables
static auto	dg::l2norm = [] ( const auto& x){ return sqrt( dg::blas1::dot(x,x));}
	Compute \( \sqrt{\sum_i x_i^2}\) using dg::blas1::dot. More...
static auto	dg::fast_l2norm
	Compute \( \sqrt{\sum_i x_i^2}\) using naive summation. More...
static auto	dg::i_control
	\( h_{n+1}= h_n \epsilon_n^{-1/p}\) More...
static auto	dg::pi_control
	\( h_{n+1}= h_n \epsilon_n^{-0.8/p}\epsilon_{n-1}^{0.31/p}\) More...
static auto	dg::pid_control
	\( h_{n+1}= h_n \epsilon_n^{-0.58/p}\epsilon_{n-1}^{0.21/p}\epsilon_{n-2}^{-0.1/p}\) More...
static auto	dg::ex_control
	\( h_{n+1} = h_n \epsilon_n^{-0.367/p}(\epsilon_n/\epsilon_{n-1})^{0.268/p} \) More...
static auto	dg::im_control
	\( h_{n+1} = h_n (h_n/h_{n-1}) \epsilon_n^{-0.98/p}(\epsilon_n/\epsilon_{n-1})^{-0.95/p} \) More...
static auto	dg::imex_control
	h_{n+1} = |ex_control| < |im_control| ? ex_control : im_control More...

Detailed Description

Enumeration Type Documentation

multistep_identifier

enum dg::multistep_identifier

Identifiers for Multistep Tableaus.

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

See also

ExplicitMultistep, ImplicitMultistep, ImExMultistep

Enumerator
IMEX_EULER_1_1
IMEX_ADAMS_2_2
IMEX_ADAMS_3_3
IMEX_KOTO_2_2
IMEX_BDF_2_2
IMEX_BDF_3_3
IMEX_BDF_4_4
IMEX_BDF_5_5
IMEX_BDF_6_6
IMEX_TVB_3_3
IMEX_TVB_4_4
IMEX_TVB_5_5
AB_1_1
AB_2_2
AB_3_3
AB_4_4
AB_5_5
eBDF_1_1
eBDF_2_2
eBDF_3_3
eBDF_4_4
eBDF_5_5
eBDF_6_6
TVB_1_1
TVB_2_2
TVB_3_3
TVB_4_4
TVB_5_5
TVB_6_6
SSP_1_1
SSP_2_2
SSP_3_2
SSP_4_2
SSP_5_3
SSP_6_3
BDF_1_1
BDF_2_2
BDF_3_3
BDF_4_4
BDF_5_5
BDF_6_6

tableau_identifier

enum dg::tableau_identifier

Identifiers for Butcher Tableaus.

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

Note

In some of the links below you might want to use the search function of your browser to find the indicated method

Enumerator
EXPLICIT_EULER_1_1	Euler
MIDPOINT_2_2	Midpoint-2-2
KUTTA_3_3	Kutta-3-3
CLASSIC_4_4	Runge-Kutta-4-4
HEUN_EULER_2_1_2	Heun-Euler-2-1-2
CAVAGLIERI_3_1_2	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems IMEXRKCB2
FEHLBERG_3_2_3	The original method uses the embedding as the solution [Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987].
FEHLBERG_4_2_3	The original method uses the embedding as the solution [Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987] (fsal)
BOGACKI_SHAMPINE_4_2_3	Bogacki-Shampine-4-2-3 (fsal)
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)
ARK324L2SA_ERK_4_2_3	ARK-4-2-3 (explicit)
ZONNEVELD_5_3_4	Zonneveld-5-3-4
ARK436L2SA_ERK_6_3_4	ARK-6-3-4 (explicit)
SAYFY_ABURUB_6_3_4	Sayfy-Aburub-6-3-4
CASH_KARP_6_4_5	Cash-Karp-6-4-5
FEHLBERG_6_4_5	Fehlberg-6-4-5
DORMAND_PRINCE_7_4_5	Dormand-Prince-7-4-5 (fsal)
TSITOURAS09_7_4_5	Tsitouras 5(4) method from 2009 (fsal), The default method in Julia
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)
ARK548L2SA_ERK_8_4_5	Kennedy and Carpenter (2019) Optimum ARK_2 method (explicit part)
VERNER_9_5_6	Verner-9-5-6 (fsal)
VERNER_10_6_7	Verner-10-6-7
FEHLBERG_13_7_8	Fehlberg-13-7-8
DORMAND_PRINCE_13_7_8	[Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987]
FEAGIN_17_8_10	Feagin-17-8-10
IMPLICIT_EULER_1_1	Euler (implicit)
IMPLICIT_MIDPOINT_1_2	implicit Midpoint
TRAPEZOIDAL_2_2	Crank-Nicolson method
SDIRK_2_1_2	generic 2nd order A and L-stable
CAVAGLIERI_IMPLICIT_3_1_2	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems
BILLINGTON_3_3_2	Billington-3-3-2
TRBDF2_3_3_2	TRBDF2-3-3-2
SANCHEZ_3_3	Sanchez-3-3
KVAERNO_4_2_3	Kvaerno-4-2-3
SDIRK_4_2_3	Cameron2002
CAVAGLIERI_IMPLICIT_4_2_3	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems The SSP scheme IMEXRKCB3c (implicit part)
ARK324L2SA_DIRK_4_2_3	ARK-4-2-3 (implicit)
SANCHEZ_3_4	Sanchez-3-4
CASH_5_2_4	Cash-5-2-4
CASH_5_3_4	Cash-5-3-4
SDIRK_5_3_4	SDIRK-5-3-4
KVAERNO_5_3_4	Kvaerno-5-3-4
ARK436L2SA_DIRK_6_3_4	ARK-6-3-4 (implicit)
KVAERNO_7_4_5	Kvaerno-7-4-5
SANCHEZ_6_5	Sanchez-6-5
ARK548L2SA_DIRK_8_4_5	Kennedy and Carpenter (2019) Optimum ARK_2 method (implicit part)
SANCHEZ_7_6	Sanchez-7-6
SSPRK_2_2	SSPRK "Shu-Osher-Form"
SSPRK_3_2	SSPRK "Shu-Osher-Form"
SSPRK_3_3	SSPRK "Shu-Osher-Form"
SSPRK_5_3	SSPRK "Shu-Osher-Form"
SSPRK_5_4	SSPRK "Shu-Osher-Form"

to

enum class dg::to

strong

Switch for the Timeloop integrate function.

Enumerator
exact	match the ending exactly
at_least	integrate to the end or further

Function Documentation

to2str()

static std::string dg::to2str ( enum to  mode )

inlinestatic

Convert integration mode to string.

Converts

• dg::to::exact to "exact"
• dg::to::at_least to "at_least"
• default "Not specified"

  Parameters

  mode

  the mode

  Returns

  string as defined above

Variable Documentation

ex_control

auto dg::ex_control

static

Initial value:

= [](auto dt, auto eps, unsigned embedded_order, unsigned order)
{
    using value_type = std::decay_t<decltype(dt[0])>;
    if( dt[1] == 0)
        return i_control( dt, eps, embedded_order, order);
    value_type m_k1 = -0.367, m_k2 = 0.268;
    value_type factor = pow( eps[0], m_k1/(value_type)embedded_order)
                      * pow( eps[0]/eps[1], m_k2/(value_type)embedded_order);
    return dt[0]*factor;
}

\( h_{n+1} = h_n \epsilon_n^{-0.367/p}(\epsilon_n/\epsilon_{n-1})^{0.268/p} \)

fast_l2norm

auto dg::fast_l2norm

static

Initial value:

= []( const auto& x){ return sqrt( dg::blas1::reduce(
            x, (double)0, dg::Sum(), dg::Square()));}

Compute \( \sqrt{\sum_i x_i^2}\) using naive summation.

The intention of this function is to be used in the Adaptive timestepping class. This function is intended for small arrays (less than 100 elements say) where the dg::blas1::dot function has its worst performance and takes overly long time to compute.

Parameters

x	Vector to take the norm of

Returns

sqrt(sum_i x_i^2) using dg::blas1::reduce

i_control

auto dg::i_control

static

Initial value:

= []( auto dt, auto eps, unsigned embedded_order, unsigned order)
{
    using value_type = std::decay_t<decltype(dt[0])>;
    return dt[0]*pow( eps[0], -1./(value_type)embedded_order);
}

\( h_{n+1}= h_n \epsilon_n^{-1/p}\)

im_control

auto dg::im_control

static

Initial value:

= []( auto dt, auto eps, unsigned embedded_order, unsigned order)
{
    using value_type = std::decay_t<decltype(dt[0])>;
    if( dt[1] == 0)
        return i_control( dt, eps, embedded_order, order);
    value_type m_k1 = -0.98, m_k2 = -0.95;
    value_type factor = pow( eps[0], m_k1/(value_type)embedded_order)
                     *  pow( eps[0]/eps[1], m_k2/(value_type)embedded_order);
    return dt[0]*dt[0]/dt[1]*factor;
}

\( h_{n+1} = h_n (h_n/h_{n-1}) \epsilon_n^{-0.98/p}(\epsilon_n/\epsilon_{n-1})^{-0.95/p} \)

imex_control

auto dg::imex_control

static

Initial value:

= []( auto dt, auto eps, unsigned embedded_order, unsigned order)
{
    using value_type = std::decay_t<decltype(dt[0])>;
    value_type h1 = ex_control( dt, eps, embedded_order, order);
    value_type h2 = im_control( dt, eps, embedded_order, order);
    
    return fabs( h1) < fabs( h2) ? h1 : h2;
}

h_{n+1} = |ex_control| < |im_control| ? ex_control : im_control

l2norm

auto dg::l2norm = [] ( const auto& x){ return sqrt( dg::blas1::dot(x,x));}

static

Compute \( \sqrt{\sum_i x_i^2}\) using dg::blas1::dot.

The intention of this function is to be used in the Adaptive timestepping class.

Parameters

x	Vector to take the norm of

Returns

sqrt(dg::blas1::dot(x,x))

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

pi_control

auto dg::pi_control

static

Initial value:

= []( auto dt, auto eps, unsigned embedded_order, unsigned order)
{
    using value_type = std::decay_t<decltype(dt[0])>;
    if( dt[1] == 0)
        return i_control( dt, eps, embedded_order, order);
    value_type m_k1 = -0.8, m_k2 = 0.31;
    value_type factor = pow( eps[0], m_k1/(value_type)embedded_order)
                     * pow( eps[1], m_k2/(value_type)embedded_order);
    return dt[0]*factor;
}

\( h_{n+1}= h_n \epsilon_n^{-0.8/p}\epsilon_{n-1}^{0.31/p}\)

pid_control

auto dg::pid_control

static

Initial value:

= []( auto dt, auto eps, unsigned embedded_order, unsigned order)
{
    using value_type = std::decay_t<decltype(dt[0])>;
    if( dt[1] == 0)
        return i_control( dt, eps, embedded_order, order);
    if( dt[2] == 0)
        return pi_control( dt, eps, embedded_order, order);
    value_type m_k1 = -0.58, m_k2 = 0.21, m_k3 = -0.1;
    
    value_type factor = pow( eps[0], m_k1/(value_type)embedded_order)
                     * pow( eps[1], m_k2/(value_type)embedded_order)
                     * pow( eps[2], m_k3/(value_type)embedded_order);
    return dt[0]*factor;
}

\( h_{n+1}= h_n \epsilon_n^{-0.58/p}\epsilon_{n-1}^{0.21/p}\epsilon_{n-2}^{-0.1/p}\)

PID stands for "Proportional" (the present error), "Integral" (the past error), "Derivative" (the future error). See a good tutorial here https://www.youtube.com/watch?v=UR0hOmjaHp0 and further information in the mathematical primer in the ARKode library. The PID controller is a good controller to start with, it does not overshoot too much, is smooth, has no systematic over- or under-estimation and converges very quickly to the desired timestep. In fact Kennedy and Carpenter, Appl. num. Math., (2003) report that it outperformed other controllers in practical problems

Template Parameters

value_type

Parameters

dt	the present (old), [0], previous [1] and second previous [2] timestep h_n
eps	the error relative to the tolerance of the present [0], previous [1] and second previous [2] timestep
embedded_order	order q of the embedded timestep
order	order p of the timestep

Returns

the new timestep