List of IMEX stabilized methods#

The PIROCK method is introduce in [AV13], the complete scheme is a IMEX scheme that allows for solving an equation of the following form (with 3 operators):

\[\dot{u} = F_R(u) + F_D(u) + F_A(u)\]

with \(F_R\) a reaction operator (solved implicitly), \(F_D\) a diffusion operator (solved explicitly by a modified ROCK2 method) and \(F_A\) an advection operator (solved with an explicit RK3 method). The first implementation of PIROCK scheme in ponio makes to solve only reaction-diffusion problem (e.g. \(F_A = 0\)), the second one to solve a complet reaction-diffusion-advection problem.

The method is divided into 5 steps:

  1. Diffusion stabilization procedure (modified ROCK2 method)

    \[\begin{split}\begin{aligned} u^{(0)} &= u^n \\ u^{(1)} &= u^n + \alpha \mu_1 \Delta t F_D(u^n) \\ u^{(j)} &= \alpha \mu_j \Delta t F_D(u^{(j-1)}) - \nu_j u^{(j-1)} - \kappa_j u^{(j-2)},\quad j=2,\dots, s-2+\ell \end{aligned}\end{split}\]

  2. Finishing procedure for diffusion

    \[\begin{split}\begin{aligned} u^{*(s-1)} &= u^{(s-2)} + \sigma_\alpha \Delta t F_D(u^{(s-2)}) \\ u^{*(s)} &= u^{*(s-1)} + \sigma_\alpha \Delta t F_D(u^{*(s-1)}) \end{aligned}\end{split}\]

  3. Starting value for advection-reaction

    \[U = u^{(s-2+\ell)}\]

  4. Finishing procedure for advection-reaction coupling

    \[\begin{split}\begin{aligned} u^{(s+1)} &= u^{(s-2+\ell)} + \gamma \Delta t F_R(u^{(s+1)}) \\ u^{(s+2)} &= u^{(s-2+\ell)} + \beta \Delta t F_D(u^{(s+1)}) + \Delta t F_A(u^{(s+1)}) + (1-2\gamma)\Delta t F_R(u^{(s+1)}) + \gamma \Delta t F_R(u^{(s+2)}) \\ u^{(s+3)} &= u^{(s-2+\ell)} + (1-2\gamma)\Delta t F_A(u^{(s+1)}) + (1-\gamma)\Delta t F_R(u^{(s+1)}) \\ u^{(s+4)} &= u^{(s-2+\ell)} + \frac{1}{3}\Delta t F_A(u^{(s+1)}) \\ u^{(s+5)} &= u^{(s-2+\ell)} + \frac{2\beta}{3} \Delta t F_D(u^{(s+1)}) + \frac{2}{3}\Delta t J_R^{-1} F_A(u^{(s+4)}) + \left(\frac{2}{3} - \gamma\right) \Delta t F_R(u^{(s+1)}) + \frac{2\gamma}{3}\Delta t F_R(u^{(s+2)}) \end{aligned}\end{split}\]

  5. Computation of the integration step

    \[\begin{split}\begin{aligned} u^{n+1} =& u^{*(s)} \\ & - \sigma_\alpha\left( 1 - \frac{\tau_\alpha}{\sigma_\alpha^2}\right)\Delta t (F_D(u^{*(s-1)}) - F_D(u^{(s-2)})) \\ & + \frac{1}{4}\Delta t F_A(u^{(s+1)}) + \frac{3}{4}\Delta t F_A(u^{(s+5)}) \\ & + \frac{1}{2}\Delta t F_R(u^{(s+1)}) + \frac{1}{2}\Delta t F_R(u^{(s+2)}) \\ & + \frac{J_R^{-\ell}}{2-4\gamma}\Delta t (F_D(u^{(s+3)}) - F_D(u^{(s+1)})) \end{aligned}\end{split}\]

    where \(\mu_j\), \(\nu_j\), \(\kappa_j\) are the same coefficients as for the standard ROCK2 method and:

    • \(\gamma = 1-\frac{\sqrt{2}}{2}\)

    • \(\beta = 1-2\alpha P'_{s-2+\ell}(0)\) where \(P'_{s-2+\ell}\) is the stability polynomial of the underlying ROCK2 method

    • \(J_R = I - \gamma \Delta t \frac{\partial F_R}{\partial u}(u^{(s-2+\ell)})\)

    • \(\sigma_\alpha = \frac{1-\alpha}{2} + \alpha \sigma\), where \(\sigma\) coming from the underlying ROCK2 method

    • \(\tau_\alpha = \frac{(\alpha-1)^2}{2} + 2\alpha(1-\alpha)\sigma + \alpha^2\tau\), where \(\sigma\) and \(\tau\) coming from the underlying ROCK2 method

For adaptive time step method, need to compute one error by operator and take the maximum value

\[\begin{split}\begin{aligned} err_D &= \sigma_\alpha \Delta t \left( 1 - \frac{\tau_\alpha}{\sigma_\alpha^2} \right) ( F_D( u^{*(s-1)} - u^{(s-2)} )) \\ err_R &= \Delta t J_R^{-1}\left( \frac{1}{6}F_R(u^{(s+1)}) - \frac{1}{6}F_R(u^{(s+2)}) \right) \\ err_A &= -\frac{3}{20}\Delta t F_A(u^{(s+1)}) + \frac{3}{10}\Delta t F_A(u^{(s+4)}) - \frac{3}{20}\Delta t F_A(u^{(s+5)}) \end{aligned}\end{split}\]

The error is computed with the maximum:

\[err = \max\left( \|err_D\|, \|err_R\|, \|err_A\|^{2/3} \right)\]

with the following error norm:

\[\|\cdot\| : err \mapsto \sum_i \left( \frac{err_i}{a_{tol} + \max(y^{n+1}_i, y^n_i) } \right)^2\]

where \(y^n\) and \(y^{n+1}\) are estimation of the solution respectively at time \(t^n\) and after an iteration \(t^{n+1}=t^n+\Delta t\), and \(a_{tol}\) and \(r_{tol}\) respectively absolute and relative tolerances. This norm give a normalized error, so its compare to \(1\) and new time step is computed as:

\[\Delta t_{new} = \sqrt{\frac{1}{err}} \Delta t\]

Note

As in a lot of adaptive time step method, the new time step stays in \(\Delta t_{new} \in [0.5\Delta t, 2 \Delta t]\) set and there is a safety factor, so the new time step is multiply by \(0.8\).

Warning

Adaptive time step method is not yet implemented for complet PIROCK ponio::runge_kutta::pirock::pirock_RDA_impl.

Still keep two free parameters \(\ell\) and \(\alpha\) given free to user, ponio provides two choices for this parameters:

  • \(\ell=2\) and \(\alpha=1\), in this case, if \(F_A=0\) and \(F_R=0\) we have the standard ROCK2 method;

  • \(\ell=1\) and \(\alpha = \frac{1}{2P'_{s-2+\ell}(0)}\), so \(\beta=0\) to minimized computation cost.

ponio provides two computers for \(\alpha\) and \(\beta\) values.

template<typename value_t = double>
class alpha_fixed#

Computer of \(\alpha\) and \(\beta\) parameters with fixed value of \(\alpha\).

Template Parameters:

value_t – Type of coefficients

Public Functions

inline value_t alpha(std::size_t, std::size_t) const#

Return \(\alpha\) (given value in constructor)

inline value_t beta(std::size_t s, std::size_t l) const#

Return \(\beta = 1 - 2\alpha P'_{s-2+\ell}(0)\).

Parameters:

  • s – number of stages

  • l – choosen parameter \(\ell\)

template<typename value_t = double>
class beta_0#

Computer of \(\alpha\) and \(\beta\) parameters with fixed value of \(\beta\) to 0.

Template Parameters:

value_t – Type of coefficients

Public Functions

inline value_t alpha(std::size_t s, std::size_t l) const#

Return \(\alpha = \frac{1}{2P'_{s-2+\ell}(0)}\) such \(\beta = 0\).

inline value_t beta(std::size_t, std::size_t) const#

Return \(\beta = 0\).

PIROCK for reaction-diffusion problem#

In this section we present the implementation of PIROCK method for only reaction-diffusion problem (i.e. \(F_A = 0\)).

template<std::size_t _l, typename alpha_beta_computer_t, typename eig_computer_t, typename _shampine_trick_caller_t = void, bool _is_embedded = false, typename _value_t = double>
class pirock_impl#

implementation of PIROCK method

Template Parameters:

  • _l – Number of augmented stages (l = 1, 2)

  • alpha_beta_computer_t – Choice of computing of \(\alpha\) and \(\beta\) parameters

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

  • _shampine_trick_caller_t – Computing of Shampine’s trick

  • _is_embedded – Set adaptive time step method (default: false)

  • value_t – Type of coefficients

Public Functions

pirock_impl() = default#

Construct a new pirock impl object.

inline pirock_impl(alpha_beta_computer_t &&_alpha_beta_computer, eig_computer_t &&_eig_computer)#

Construct a new pirock impl object.

Parameters:

  • _alpha_beta_computer – computer of parameters alpha and beta object that has two member functions which take number of stages (s) l parameter

  • _eig_computer – eigenvalue computer functor (that take as argument the function of explicit part, the current time, the current state and the current time step)

template<typename _shampine_trick_caller_t_>
inline pirock_impl(alpha_beta_computer_t &&_alpha_beta_computer, eig_computer_t &&_eig_computer, _shampine_trick_caller_t_ &&_shampine_trick_caller)#

Construct a new pirock impl object with Shampine’s trick.

Parameters:

  • _alpha_beta_computer – alpha and beta computer object

  • _eig_computer – eigenvalue computer functor

  • _shampine_trick_caller – Shampine’s trick functor

template<typename problem_t, typename state_t, typename array_ki_t>
inline void operator()(problem_t &pb, value_t &tn, state_t &un, array_ki_t &U, value_t &dt, state_t &u_np1)#

iteration of PIROCK method

Template Parameters:

  • problem_t – type of \(f\)

  • state_t – type of current state

  • array_ki_t – type of temporary stages (only 3 needed for ROCK2)

Parameters:

  • pb – problem \((F_D, F_R)\) and the Jacibian of reaction part \(\frac{\partial F_R}{\partial u}\)

  • tn – current time

  • un – current state

  • U – array of temporary stages

  • dt – current time step

  • u_np1 – solution \(u^{n+1}\) at time \(t^{n+1} = t^n + \Delta t\)

inline auto &info()#

gets iteration_info object

inline auto const &info() const#

gets iteration_info object (constant version)

template<typename rock_t = rock_coeff>
inline auto &abs_tol(value_t tol_)#

set absolute tolerance in chained config

Parameters:

tol_ – tolerance

Returns:

auto& returns this object

template<typename rock_t = rock_coeff>
inline auto &rel_tol(value_t tol_)#

set relative tolerance in chained config

Parameters:

tol_ – tolerance

Returns:

auto& returns this object

Helper functions#

User interface functions to build a PIROCK method.

template<std::size_t l = 1, bool is_embedded = false, typename value_t = double, typename alpha_beta_computer_t, typename eig_computer_t, typename shampine_trick_caller_t>
auto ponio::runge_kutta::pirock::pirock(alpha_beta_computer_t &&alpha_beta_computer, eig_computer_t &&eig_computer, shampine_trick_caller_t &&shampine_trick_caller)#

helper function to build PIROCK algorithm

Template Parameters:

  • l – Number of augmented stages (l = 1, 2)

  • is_embedded – Set adaptive time step method (default: false)

  • value_t – Type of coefficients

  • alpha_beta_computer_t – Choice of computing of \(\alpha\) and \(\beta\) parameters

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

  • shampine_trick_caller_t – Computing of Shampine’s trick

Parameters:

  • alpha_beta_computer – \(\alpha\) and \(\beta\) computer object

  • eig_computer – Eigenvalue computer of explicit part of the problem

  • shampine_trick_caller – Shampine’s trick computer

template<std::size_t l = 1, typename value_t = double, typename alpha_beta_computer_t, typename eig_computer_t>
auto ponio::runge_kutta::pirock::pirock(alpha_beta_computer_t &&alpha_beta_computer, eig_computer_t &&eig_computer)#

helper function to build PIROCK algorithm

Note

Without a Shampine’s trick caller, this method is fixed time step.

Template Parameters:

  • l – Number of augmented stages (l = 1, 2)

  • value_t – Type of coefficients

  • alpha_beta_computer_t – Choice of computing of \(\alpha\) and \(\beta\) parameters

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

Parameters:

  • alpha_beta_computer – \(\alpha\) and \(\beta\) computer object

  • eig_computer – Eigenvalue computer of explicit part of the problem

template<std::size_t l = 1, typename value_t = double, typename eig_computer_t>
auto ponio::runge_kutta::pirock::pirock(eig_computer_t &&eig_computer)#

helper function to build PIROCK algorithm

Note

Without a Shampine’s trick caller, this method is fixed time step, and without a \(\alpha\) and \(\beta\) computer, parameters are fixed to \(\beta = 0\).

Template Parameters:

  • l – Number of augmented stages (l = 1, 2)

  • value_t – Type of coefficients

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

Parameters:

eig_computer – Eigenvalue computer of explicit part of the problem

template<std::size_t l = 1, typename value_t = double>
auto ponio::runge_kutta::pirock::pirock()#

helper function to build PIROCK algorithm

Note

Without a Shampine’s trick caller, this method is fixed time step, without a \(\alpha\) and \(\beta\) computer, parameters are fixed to \(\beta = 0\) and without a eigenvalues computer this method use power method to estimate spectral radius of explicit part (diffusion operator).

Template Parameters:

  • l – Number of augmented stages (l = 1, 2)

  • value_t – Type of coefficients

\(\ell=2\) and \(\alpha = 1\) case#

Following functions are useful for to build a PIROCK method with \(\ell=2\) and \(\alpha = 1\) (with ponio::runge_kutta::pirock::alpha_fixed computer).

template<typename value_t = double, typename eig_computer_t>
auto ponio::runge_kutta::pirock::pirock_a1(eig_computer_t &&eig_computer)#

helper function to build PIROCK algorithm

Note

This function builds a PIROCK algorithm with parameters \(\ell=2\) and \(\alpha = 1\) and specific spectral radius estimator of explicit part (diffusion operator).

Template Parameters:

  • value_t – Type of coefficients

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

Parameters:

eig_computer – Eigenvalue computer of explicit part of the problem

template<typename value_t = double>
auto ponio::runge_kutta::pirock::pirock_a1()#

helper function to build PIROCK algorithm

Note

This function builds a PIROCK algorithm with parameters \(\ell=2\) and \(\alpha = 1\) and power method to estimate spectral radius.

Template Parameters:

value_t – Type of coefficients

\(\ell=1\) and \(\beta = 0\) case#

Following functions are useful for to build a PIROCK method with \(\ell=1\) and \(\beta = 0\) (with ponio::runge_kutta::pirock::beta_0 computer).

template<typename value_t = double, typename eig_computer_t>
auto ponio::runge_kutta::pirock::pirock_b0(eig_computer_t &&eig_computer)#

helper function to build PIROCK algorithm

Note

This function builds a PIROCK algorithm with parameters \(\ell=1\) and \(\beta = 0\) and specific spectral radius estimator of explicit part (diffusion operator).

Template Parameters:

  • value_t – Type of coefficients

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

Parameters:

eig_computer – Eigenvalue computer of explicit part of the problem

template<typename value_t = double>
auto ponio::runge_kutta::pirock::pirock_b0()#

helper function to build PIROCK algorithm

Note

This function builds a PIROCK algorithm with parameters \(\ell=1\) and \(\beta = 0\) and power method to estimate spectral radius.

Template Parameters:

value_t – Type of coefficients

PIROCK for reaction-diffusion-advection problem#

In this section we present the implementation of complet PIROCK method (i.e. for reaction-diffusion-advection problem).

template<std::size_t _l, typename alpha_beta_computer_t, typename eig_computer_t, typename _shampine_trick_caller_t = void, bool _is_embedded = false, typename _value_t = double>
class pirock_RDA_impl#

implementation of PIROCK method for reaction-diffusion-advection problem

Template Parameters:

  • _l – Number of augmented stages (l = 1, 2)

  • alpha_beta_computer_t – Choice of computing of \(\alpha\) and \(\beta\) parameters

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

  • _shampine_trick_caller_t – Computing of Shampine’s trick

  • _is_embedded – Set adaptive time step method (default: false)

  • value_t – Type of coefficients

Public Functions

pirock_RDA_impl() = default#

Construct a new pirock_RDA_impl object.

inline pirock_RDA_impl(alpha_beta_computer_t &&_alpha_beta_computer, eig_computer_t &&_eig_computer)#

Construct a new pirock_RDA_impl object.

Parameters:

  • _alpha_beta_computer – computer of parameters alpha and beta object that has two member functions which take number of stages (s) l parameter

  • _eig_computer – eigenvalue computer functor (that take as argument the function of explicit part, the current time, the current state and the current time step)

template<typename _shampine_trick_caller_t_>
inline pirock_RDA_impl(alpha_beta_computer_t &&_alpha_beta_computer, eig_computer_t &&_eig_computer, _shampine_trick_caller_t_ &&_shampine_trick_caller)#

Construct a new pirock_RDA_impl object with Shampine’s trick.

Parameters:

  • _alpha_beta_computer – alpha and beta computer object

  • _eig_computer – eigenvalue computer functor

  • _shampine_trick_caller – Shampine’s trick functor

template<typename problem_t, typename state_t, typename array_ki_t>
inline void operator()(problem_t &pb, value_t &tn, state_t &un, array_ki_t &U, value_t &dt, state_t &u_np1)#

iteration of PIROCK_RDA method

Template Parameters:

  • problem_t – type of \(f\)

  • state_t – type of current state

  • array_ki_t – type of temporary stages (13 or 18 stages needed)

Parameters:

  • pb – problem with 3 operators: \(F_R\), \(F_D\) and \(F_A\)

  • tn – current time

  • un – current state

  • U – array of temporary stages

  • dt – current time step

  • u_np1 – solution \(u^{n+1}\) at time \(t^{n+1} = t^n + \Delta t\)

inline auto &info()#

gets iteration_info object

inline auto const &info() const#

gets iteration_info object (constant version)

template<typename rock_t = rock_coeff>
inline auto &abs_tol(value_t tol_)#

set absolute tolerance in chained config

Parameters:

tol_ – tolerance

Returns:

auto& returns this object

template<typename rock_t = rock_coeff>
inline auto &rel_tol(value_t tol_)#

set relative tolerance in chained config

Parameters:

tol_ – tolerance

Returns:

auto& returns this object

Helper functions#

User interface functions to build a PIROCK method.

template<std::size_t l = 1, bool is_embedded = false, typename value_t = double, typename alpha_beta_computer_t, typename eig_computer_t, typename shampine_trick_caller_t>
auto ponio::runge_kutta::pirock::pirock_RDA(alpha_beta_computer_t &&alpha_beta_computer, eig_computer_t &&eig_computer, shampine_trick_caller_t &&shampine_trick_caller)#

helper function to build PIROCK algorithm

Template Parameters:

  • l – Number of augmented stages (l = 1, 2)

  • is_embedded – Set adaptive time step method (default: false)

  • value_t – Type of coefficients

  • alpha_beta_computer_t – Choice of computing of \(\alpha\) and \(\beta\) parameters

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

  • shampine_trick_caller_t – Computing of Shampine’s trick

Parameters:

  • alpha_beta_computer – \(\alpha\) and \(\beta\) computer object

  • eig_computer – Eigenvalue computer of explicit part of the problem

  • shampine_trick_caller – Shampine’s trick computer

template<std::size_t l = 1, typename value_t = double, typename alpha_beta_computer_t, typename eig_computer_t>
auto ponio::runge_kutta::pirock::pirock_RDA(alpha_beta_computer_t &&alpha_beta_computer, eig_computer_t &&eig_computer)#

helper function to build PIROCK algorithm

Note

Without a Shampine’s trick caller, this method is fixed time step.

Template Parameters:

  • l – Number of augmented stages (l = 1, 2)

  • value_t – Type of coefficients

  • alpha_beta_computer_t – Choice of computing of \(\alpha\) and \(\beta\) parameters

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

Parameters:

  • alpha_beta_computer – \(\alpha\) and \(\beta\) computer object

  • eig_computer – Eigenvalue computer of explicit part of the problem

template<std::size_t l = 1, typename value_t = double, typename eig_computer_t>
auto ponio::runge_kutta::pirock::pirock_RDA(eig_computer_t &&eig_computer)#

helper function to build PIROCK algorithm

Note

Without a Shampine’s trick caller, this method is fixed time step, and without a \(\alpha\) and \(\beta\) computer, parameters are fixed to \(\beta = 0\).

Template Parameters:

  • l – Number of augmented stages (l = 1, 2)

  • value_t – Type of coefficients

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

Parameters:

eig_computer – Eigenvalue computer of explicit part of the problem

template<std::size_t l = 1, typename value_t = double>
auto ponio::runge_kutta::pirock::pirock_RDA()#

helper function to build PIROCK algorithm

Note

Without a Shampine’s trick caller, this method is fixed time step, without a \(\alpha\) and \(\beta\) computer, parameters are fixed to \(\beta = 0\) and without a eigenvalues computer this method use power method to estimate spectral radius of explicit part (diffusion operator).

Template Parameters:

  • l – Number of augmented stages (l = 1, 2)

  • value_t – Type of coefficients

\(\ell=2\) and \(\alpha = 1\) case#

Following functions are useful for to build a PIROCK method with \(\ell=2\) and \(\alpha = 1\) (with ponio::runge_kutta::pirock::alpha_fixed computer).

template<typename value_t = double, typename eig_computer_t>
auto ponio::runge_kutta::pirock::pirock_RDA_a1(eig_computer_t &&eig_computer)#

helper function to build PIROCK algorithm

Note

This function builds a PIROCK algorithm with parameters \(\ell=2\) and \(\alpha = 1\) and specific spectral radius estimator of explicit part (diffusion operator).

Template Parameters:

  • value_t – Type of coefficients

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

Parameters:

eig_computer – Eigenvalue computer of explicit part of the problem

template<typename value_t = double>
auto ponio::runge_kutta::pirock::pirock_RDA_a1()#

helper function to build PIROCK algorithm

Note

This function builds a PIROCK algorithm with parameters \(\ell=2\) and \(\alpha = 1\) and power method to estimate spectral radius.

Template Parameters:

value_t – Type of coefficients

\(\ell=1\) and \(\beta = 0\) case#

Following functions are useful for to build a PIROCK method with \(\ell=1\) and \(\beta = 0\) (with ponio::runge_kutta::pirock::beta_0 computer).

template<typename value_t = double, typename eig_computer_t>
auto ponio::runge_kutta::pirock::pirock_RDA_b0(eig_computer_t &&eig_computer)#

helper function to build PIROCK algorithm

Note

This function builds a PIROCK algorithm with parameters \(\ell=1\) and \(\beta = 0\) and specific spectral radius estimator of explicit part (diffusion operator).

Template Parameters:

  • value_t – Type of coefficients

  • eig_computer_t – Computing of eigenvalues of explicit part (diffusion)

Parameters:

eig_computer – Eigenvalue computer of explicit part of the problem

template<typename value_t = double>
auto ponio::runge_kutta::pirock::pirock_RDA_b0()#

helper function to build PIROCK algorithm

Note

This function builds a PIROCK algorithm with parameters \(\ell=1\) and \(\beta = 0\) and power method to estimate spectral radius.

Template Parameters:

value_t – Type of coefficients