Reference documentation for deal.II version 9.2.0
SUNDIALS::ARKode< VectorType > Class Template Reference

#include <deal.II/sundials/arkode.h>

## Public Member Functions

~ARKode ()

unsigned int solve_ode (VectorType &solution)

void reset (const double t, const double h, const VectorType &y)

DeclException1 (ExcARKodeError, int,<< "One of the SUNDIALS ARKode internal functions "<< " returned a negative error code: "<< arg1<< ". Please consult SUNDIALS manual.")

## Public Attributes

std::function< void(VectorType &)> reinit_vector

std::function< int(const double t, const VectorType &y, VectorType &explicit_f)> explicit_function

std::function< int(const double t, const VectorType &y, VectorType &res)> implicit_function

std::function< int(const int convfail, const double t, const double gamma, const VectorType &ypred, const VectorType &fpred, bool &j_is_current)> setup_jacobian

std::function< int(const double t, const double gamma, const VectorType &ycur, const VectorType &fcur, const VectorType &rhs, VectorType &dst)> solve_jacobian_system

std::function< int(const double t)> setup_mass

std::function< int(const VectorType &rhs, VectorType &dst)> solve_mass_system

std::function< void(const double t, const VectorType &sol, const unsigned int step_number)> output_step

std::function< bool(const double t, VectorType &sol)> solver_should_restart

std::function< VectorType &()> get_local_tolerances

## Private Member Functions

DeclException1 (ExcFunctionNotProvided, std::string,<< "Please provide an implementation for the function << arg1<< "\")

void set_functions_to_trigger_an_assert ()

## Private Attributes

void * arkode_mem

N_Vector yy

N_Vector abs_tolls

MPI_Comm communicator

GrowingVectorMemory< VectorTypemem

## Detailed Description

### template<typename VectorType = Vector<double>> class SUNDIALS::ARKode< VectorType >

Interface to SUNDIALS additive Runge-Kutta methods (ARKode).

The class ARKode is a wrapper to SUNDIALS variable-step, embedded, additive Runge-Kutta solver which is a general purpose solver for systems of ordinary differential equations characterized by the presence of both fast and slow dynamics.

Fast dynamics are treated implicitly, and slow dynamics are treated explicitly, using nested families of implicit and explicit Runge-Kutta solvers.

Citing directly from ARKode documentation:

ARKode solves ODE initial value problems (IVPs) in $$R^N$$. These problems should be posed in explicit form as

$M\dot y = f_E(t, y) + f_I (t, y), \qquad y(t_0) = y_0.$

Here, $$t$$ is the independent variable (e.g. time), and the dependent variables are given by $$y \in R^N$$, and we use notation $$\dot y$$ to denote $$dy/dt$$. $$M$$ is a user-supplied nonsingular operator from $$R^N \to R^N$$. This operator may depend on $$t$$ but not on $$y$$.

For standard systems of ordinary differential equations and for problems arising from the spatial semi-discretization of partial differential equations using finite difference or finite volume methods, $$M$$ is typically the identity matrix, $$I$$. For PDEs using a finite-element spatial semi-discretization $$M$$ is typically a well-conditioned mass matrix.

The two right-hand side functions may be described as:

• $$f_E(t, y)$$: contains the "slow" time scale components of the system. This will be integrated using explicit methods.
• $$f_I(t, y)$$: contains the "fast" time scale components of the system. This will be integrated using implicit methods.

ARKode may be used to solve stiff, nonstiff and multi-rate problems. Roughly speaking, stiffness is characterized by the presence of at least one rapidly damped mode, whose time constant is small compared to the time scale of the solution itself. In the implicit/explicit (ImEx) splitting above, these stiff components should be included in the right-hand side function $$f_I (t, y)$$.

For multi-rate problems, a user should provide both of the functions $$f_E$$ and $$f_I$$ that define the IVP system.

For nonstiff problems, only $$f_E$$ should be provided, and $$f_I$$ is assumed to be zero, i.e. the system reduces to the non-split IVP:

$M\dot y = f_E(t, y), \qquad y(t_0) = y_0.$

In this scenario, the ARK methods reduce to classical explicit Runge-Kutta methods (ERK). For these classes of methods, ARKode allows orders of accuracy $$q = \{2, 3, 4, 5, 6, 8\}$$, with embeddings of orders $$p = \{1, 2, 3, 4, 5, 7\}$$. These default to the Heun-Euler-2-1-2, Bogacki-Shampine-4-2-3, Zonneveld-5-3-4, Cash-Karp-6-4-5, Verner-8-5-6 and Fehlberg-13-7-8 methods, respectively.

Finally, for stiff (linear or nonlinear) problems the user may provide only $$f_I$$, implying that $$f_E = 0$$, so that the system reduces to the non-split IVP

$M\dot y = f_I(t, y), \qquad y(t_0) = y_0.$

Similarly to ERK methods, in this scenario the ARK methods reduce to classical diagonally-implicit Runge-Kutta methods (DIRK). For these classes of methods, ARKode allows orders of accuracy $$q = \{2, 3, 4, 5\}$$, with embeddings of orders $$p = \{1, 2, 3, 4\}$$. These default to the SDIRK-2-1-2, ARK-4-2-3 (implicit), SDIRK-5-3-4 and ARK-8-4-5 (implicit) methods, respectively.

For both DIRK and ARK methods, an implicit system of the form

$G(z_i) \dealcoloneq M z_i - h_n A^I_{i,i} f_I (t^I_{n,i}, z_i) - a_i = 0$

must be solved for each stage $$z_i , i = 1, \ldots, s$$, where we have the data

$a_i \dealcoloneq M y_{n-1} + h_n \sum_{j=1}^{i-1} [ A^E_{i,j} f_E(t^E_{n,j}, z_j) + A^I_{i,j} f_I (t^I_{n,j}, z_j)]$

for the ARK methods, or

$a_i \dealcoloneq M y_{n-1} + h_n \sum_{j=1}^{i-1} A^I_{i,j} f_I (t^I_{n,j}, z_j)$

for the DIRK methods. Here $$A^I_{i,j}$$ and $$A^E_{i,j}$$ are the Butcher's tables for the chosen solver.

If $$f_I(t,y)$$ depends nonlinearly on $$y$$ then the systems above correspond to a nonlinear system of equations; if $$f_I (t, y)$$ depends linearly on $$y$$ then this is a linear system of equations. By specifying the flag implicit_function_is_linear, ARKode takes some shortcuts that allow a faster solution process.

For systems of either type, ARKode allows a choice of solution strategy. The default solver choice is a variant of Newton's method,

$z_i^{m+1} = z_i^m +\delta^{m+1},$

where $$m$$ is the Newton index, and the Newton update $$\delta^{m+1}$$ requires the solution of the linear Newton system

$N(z_i^m) \delta^{m+1} = -G(z_i^m),$

where

$N \dealcoloneq M - \gamma J, \quad J \dealcoloneq \frac{\partial f_I}{\partial y}, \qquad \gamma\dealcoloneq h_n A^I_{i,i}.$

As an alternate to Newton's method, ARKode may solve for each stage $$z_i ,i = 1, \ldots , s$$ using an Anderson-accelerated fixed point iteration

$z_i^{m+1} = g(z_i^{m}), m=0,1,\ldots.$

Unlike with Newton's method, this option does not require the solution of a linear system at each iteration, instead opting for solution of a low-dimensional least-squares solution to construct the nonlinear update.

Finally, if the user specifies implicit_function_is_linear, i.e., $$f_I(t, y)$$ depends linearly on $$y$$, and if the Newton-based nonlinear solver is chosen, then the system will be solved using only a single Newton iteration. Notice that in order for the Newton solver to be used, at least the solve_jacobian_system() function should be supplied. If this function is not supplied, then only the fixed-point iteration will be supported, and the implicit_function_is_linear setting is ignored.

The optimal solver (Newton vs fixed-point) is highly problem-dependent. Since fixed-point solvers do not require the solution of any linear systems, each iteration may be significantly less costly than their Newton counterparts. However, this can come at the cost of slower convergence (or even divergence) in comparison with Newton-like methods. These fixed-point solvers do allow for user specification of the Anderson-accelerated subspace size, $$m_k$$. While the required amount of solver memory grows proportionately to $$m_k N$$, larger values of $$m_k$$ may result in faster convergence.

This improvement may be significant even for "small" values, e.g. $$1 \leq m_k \leq 5$$, and convergence may not improve (or even deteriorate) for larger values of $$m_k$$. While ARKode uses a Newton-based iteration as its default solver due to its increased robustness on very stiff problems, it is highly recommended that users also consider the fixed-point solver for their cases when attempting a new problem.

For either the Newton or fixed-point solvers, it is well-known that both the efficiency and robustness of the algorithm intimately depends on the choice of a good initial guess. In ARKode, the initial guess for either nonlinear solution method is a predicted value $$z_i(0)$$ that is computed explicitly from the previously-computed data (e.g. $$y_{n-2}, y_{n-1}$$, and $$z_j$$ where $$j < i$$). Additional information on the specific predictor algorithms implemented in ARKode is provided in ARKode documentation.

The user has to provide the implementation of the following std::function:

• reinit_vector; and either one or both of
• implicit_function;
• explicit_function;

If the mass matrix is different from the identity, the user should supply

• solve_mass_system; and, optionally,
• setup_mass;

If the use of a Newton method is desired, then the user should also supply

• solve_jacobian_system; and optionally
• setup_jacobian;

Also the following functions could be rewritten. By default they do nothing, or are not required.

• solver_should_restart;
• get_local_tolerances;

To produce output at fixed steps, overload the function

• output_step;

To provide a simple example, consider the harmonic oscillator problem:

$\begin{split} u'' & = -k^2 u \\ u (0) & = 0 \\ u'(0) & = k \end{split}$

We write it in terms of a first order ode:

$\begin{matrix} y_0' & = y_1 \\ y_1' & = - k^2 y_0 \end{matrix}$

That is $$y' = A y$$ where

$A \dealcoloneq \begin{matrix} 0 & 1 \\ -k^2 &0 \end{matrix}$

and $$y(0)=(0, k)$$.

The exact solution is $$y_0(t) = \sin(k t)$$, $$y_1(t) = y_0'(t) = k \cos(k *t)$$, $$y_1'(t) = -k^2 \sin(k t)$$.

A minimal implementation, using only explicit RK methods, is given by the following code snippet:

{
v.reinit(2);
};
double kappa = 1.0;
ode.explicit_function = [kappa] (double,
const VectorType &y,
VectorType &ydot) -> int
{
ydot[0] = y[1];
ydot[1] = -kappa*kappa*y[0];
return 0;
};
y[1] = kappa;
ode.solve_ode(y);

Definition at line 313 of file arkode.h.

## ◆ ARKode()

template<typename VectorType >
 SUNDIALS::ARKode< VectorType >::ARKode ( const AdditionalData & data = AdditionalData(), const MPI_Comm mpi_comm = MPI_COMM_WORLD )

Constructor. It is possible to fine tune the SUNDIALS ARKode solver by passing an AdditionalData() object that sets all of the solver parameters.

The MPI communicator is simply ignored in the serial case.

Parameters
 data ARKode configuration data mpi_comm MPI communicator

Definition at line 239 of file arkode.cc.

## ◆ ~ARKode()

template<typename VectorType >
 SUNDIALS::ARKode< VectorType >::~ARKode ( )

Destructor.

Definition at line 253 of file arkode.cc.

## ◆ solve_ode()

template<typename VectorType >
 unsigned int SUNDIALS::ARKode< VectorType >::solve_ode ( VectorType & solution )

Integrate the initial value problem. This function returns the final number of computed steps.

Definition at line 271 of file arkode.cc.

## ◆ reset()

template<typename VectorType >
 void SUNDIALS::ARKode< VectorType >::reset ( const double t, const double h, const VectorType & y )

Clear internal memory and start with clean objects. This function is called when the simulation starts and when the user returns true to a call to solver_should_restart().

By default solver_should_restart() returns false. If the user needs to implement, for example, local adaptivity in space, he or she may assign a different function to solver_should_restart() that performs all mesh changes, transfers the solution to the new mesh, and returns true.

Parameters
 [in] t The new starting time [in] h The new starting time step [in,out] y The new initial solution

Definition at line 353 of file arkode.cc.

## ◆ DeclException1() [1/2]

template<typename VectorType = Vector<double>>
 SUNDIALS::ARKode< VectorType >::DeclException1 ( ExcARKodeError , int , << "One of the SUNDIALS ARKode< VectorType > internal functions "<< " returned a negative error code: "<< arg1<< ". Please consult SUNDIALS manual." )

Handle ARKode exceptions.

## ◆ DeclException1() [2/2]

template<typename VectorType = Vector<double>>
 SUNDIALS::ARKode< VectorType >::DeclException1 ( ExcFunctionNotProvided , std::string , << "Please provide an implementation for the function \ arg1<< "\" )
private

Throw an exception when a function with the given name is not implemented.

## ◆ set_functions_to_trigger_an_assert()

template<typename VectorType >
 void SUNDIALS::ARKode< VectorType >::set_functions_to_trigger_an_assert ( )
private

This function is executed at construction time to set the std::function above to trigger an assert if they are not implemented.

Definition at line 494 of file arkode.cc.

## ◆ reinit_vector

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::reinit_vector

A function object that users need to supply and that is intended to reinit the given vector.

Definition at line 549 of file arkode.h.

## ◆ explicit_function

template<typename VectorType = Vector<double>>
 std::function< int(const double t, const VectorType &y, VectorType &explicit_f)> SUNDIALS::ARKode< VectorType >::explicit_function

A function object that users may supply and that is intended to compute the explicit part of the IVP right hand side. Sets $$explicit_f = f_E(t, y)$$.

At least one of explicit_function() or implicit_function() must be provided. According to which one is provided, explicit, implicit, or mixed RK methods are used.

This function should return:

• 0: Success
• >0: Recoverable error (ARKodeReinit will be called if this happens, and then last function will be attempted again
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 569 of file arkode.h.

## ◆ implicit_function

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::implicit_function

A function object that users may supply and that is intended to compute the implicit part of the IVP right hand side. Sets $$implicit_f = f_I(t, y)$$.

At least one of explicit_function() or implicit_function() must be provided. According to which one is provided, explicit, implicit, or mixed RK methods are used.

This function should return:

• 0: Success
• >0: Recoverable error (ARKodeReinit will be called if this happens, and then last function will be attempted again
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 588 of file arkode.h.

## ◆ setup_jacobian

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::setup_jacobian

A function object that users may supply and that is intended to prepare the linear solver for subsequent calls to solve_jacobian_system().

Make sure that after a call to this function, we know how to compute solutions of systems $$A x = b$$, where $$A$$ is some approximation to the Newton matrix, $$M - \gamma \partial f_I/\partial y$$. This function is optional. If the user does not provide it, then solve_jacobian_system() is assumed to also perform the setup internally.

The setup_jacobian() function may call a user-supplied function to compute needed data related to the Jacobian matrix. Alternatively, it may choose to retrieve and use stored values of this data. In either case, setup_jacobian() may also preprocess that data as needed for solve_jacobian_system(), which may involve calling a generic function (such as for LU factorization).

This data may be intended either for direct use (in a direct linear solver) or for use in a preconditioner (in a preconditioned iterative linear solver). The setup_jacobian() function is not called at every stage solve (or even every time step), but only as frequently as the solver determines that it is appropriate to perform the setup task. In this way, Jacobian-related data generated by setup_jacobian() is expected to be used over a number of time steps.

If the user uses a matrix based computation of the Jacobian, then this is the right place where an assembly routine shoulde be called to assemble both a matrix and a preconditioner for the Jacobian system. Subsequent calls (possibly more than one) to solve_jacobian_system() can assume that this function has been called at least once.

Notice that no assumption is made by this interface on what the user should do in this function. ARKode only assumes that after a call to setup_jacobian() it is possible to call solve_jacobian_system(), to obtain a solution $$x$$ to the system $$J x = b$$. If this function is not provided, then it is never called.

Arguments to the function are

Parameters
 [in] t the current time [in] gamma the current factor to use in the jacobian computation [in] ypred is the predicted $$y$$ vector for the current ARKode internal step [in] fpred is the value of the implicit right-hand side at ypred, $$f_I (t_n, ypred)$$. [in] convfail Input flag used to indicate any problem that occurred during the solution of the nonlinear equation on the current time step for which the linear solver is being used. This flag can be used to help decide whether the Jacobian data kept by a linear solver needs to be updated or not. Its possible values are:
• ARK_NO_FAILURES: this value is passed if either this is the first call for this step, or the local error test failed on the previous attempt at this step (but the Newton iteration converged).
• ARK_FAIL_BAD_J: this value is passed if (a) the previous Newton corrector iteration did not converge and the linear solver's setup function indicated that its Jacobian-related data is not current, or (b) during the previous Newton corrector iteration, the linear solver's solve function failed in a recoverable manner and the linear solver's setup function indicated that its Jacobian-related data is not current.
• ARK_FAIL_OTHER: this value is passed if during the current internal step try, the previous Newton iteration failed to converge even though the linear solver was using current Jacobian-related data.
Parameters
 [out] j_is_current a boolean to be filled in by setup_jacobian(). The value should be set to true if the Jacobian data is current after the call, and should be set to false if its Jacobian data is not current. If setup_jacobian() calls for re-evaluation of Jacobian data (based on convfail and ARKode state data), then it should set j_is_current to true unconditionally, otherwise an infinite loop can result.

This function should return:

• 0: Success
• >0: Recoverable error (ARKodeReinit will be called if this happens, and then last function will be attempted again
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 680 of file arkode.h.

## ◆ solve_jacobian_system

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::solve_jacobian_system

A function object that users may supply and that is intended to solve the Jacobian linear system. This function will be called by ARKode (possibly several times) after setup_jacobian() has been called at least once. ARKode tries to do its best to call setup_jacobian() the minimum amount of times. If convergence can be achieved without updating the Jacobian, then ARKode does not call setup_jacobian() again. If, on the contrary, internal ARKode convergence tests fail, then ARKode calls again setup_jacobian() with updated vectors and coefficients so that successive calls to solve_jacobian_systems() lead to better convergence in the Newton process.

If you do not specify a solve_jacobian_system() function, then a fixed point iteration is used instead of a Newton method. Notice that this may not converge, or may converge very slowly.

The jacobian $$J$$ should be (an approximation of) the system Jacobian

$J = M - \gamma \frac{\partial f_I}{\partial y}$

evaluated at t, ycur. fcur is $$f_I(t,ycur)$$.

A call to this function should store in dst the result of $$J^{-1}$$ applied to src, i.e., J*dst = src. It is the users responsibility to set up proper solvers and preconditioners inside this function.

Arguments to the function are

Parameters
 [in] t the current time [in] gamma the current factor to use in the jacobian computation [in] ycur is the current $$y$$ vector for the current ARKode internal step [in] fcur is the current value of the implicit right-hand side at ycur, $$f_I (t_n, ypred)$$.

This function should return:

• 0: Success
• >0: Recoverable error (ARKodeReinit will be called if this happens, and then last function will be attempted again
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 732 of file arkode.h.

## ◆ setup_mass

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::setup_mass

A function object that users may supply and that is intended to setup the mass matrix. This function is called by ARKode any time a mass matrix update is required. The user should compute the mass matrix (or update all the variables that allow the application of the mass matrix). This function is called by ARKode once, before any call to solve_mass_system().

ARKode supports the case where the mass matrix may depend on time, but not the case where the mass matrix depends on the solution itself.

If the user does not provide a solve_mass_matrix() function, then the identity is used. If the setup_mass() function is not provided, then solve_mass_system() should do all the work by itself.

If the user uses a matrix based computation of the mass matrix, then this is the right place where an assembly routine shoulde be called to assemble both a matrix and a preconditioner. Subsequent calls (possibly more than one) to solve_mass_system() can assume that this function has been called at least once.

Notice that no assumption is made by this interface on what the user should do in this function. ARKode only assumes that after a call to setup_mass() it is possible to call solve_mass_system(), to obtain a solution $$x$$ to the system $$M x = b$$.

This function should return:

• 0: Success
• >0: Recoverable error (ARKodeReinit will be called if this happens, and then last function will be attempted again
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 768 of file arkode.h.

## ◆ solve_mass_system

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::solve_mass_system

A function object that users may supply and that is intended to solve the mass matrix linear system. This function will be called by ARKode (possibly several times) after setup_mass() has been called at least once. ARKode tries to do its best to call setup_mass() the minimum amount of times.

A call to this function should store in dst the result of $$M^{-1}$$ applied to src, i.e., M*dst = src. It is the users responsibility to set up proper solvers and preconditioners inside this function.

This function should return:

• 0: Success
• >0: Recoverable error (ARKodeReinit will be called if this happens, and then last function will be attempted again
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 789 of file arkode.h.

## ◆ output_step

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::output_step

A function object that users may supply and that is intended to postprocess the solution. This function is called by ARKode at fixed time increments (every output_period seconds), and it is passed a polynomial interpolation of the solution, computed using the current ARK order and the (internally stored) previously computed solution steps.

Notice that it is well possible that internally ARKode computes a time step which is much larger than the output_period step, and therefore calls this function consecutively several times by simply performing all intermediate interpolations. There is no relationship between how many times this function is called and how many time steps have actually been computed.

Definition at line 808 of file arkode.h.

## ◆ solver_should_restart

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::solver_should_restart

A function object that users may supply and that is intended to evaluate whether the solver should be restarted (for example because the number of degrees of freedom has changed).

This function is supposed to perform all operations that are necessary in sol to make sure that the resulting vectors are consistent, and of the correct final size.

For example, one may decide that a local refinement is necessary at time t. This function should then return true, and change the dimension of sol to reflect the new dimension. Since ARKode does not know about the new dimension, an internal reset is necessary.

The default implementation simply returns false, i.e., no restart is performed during the evolution.

Definition at line 827 of file arkode.h.

## ◆ get_local_tolerances

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::ARKode< VectorType >::get_local_tolerances

A function object that users may supply and that is intended to return a vector whose components are the weights used by ARKode to compute the vector norm. The implementation of this function is optional, and it is used only if implemented.

Definition at line 835 of file arkode.h.

## ◆ data

template<typename VectorType = Vector<double>>
private

ARKode configuration data.

Definition at line 868 of file arkode.h.

## ◆ arkode_mem

template<typename VectorType = Vector<double>>
 void* SUNDIALS::ARKode< VectorType >::arkode_mem
private

ARKode memory object.

Definition at line 873 of file arkode.h.

## ◆ yy

template<typename VectorType = Vector<double>>
 N_Vector SUNDIALS::ARKode< VectorType >::yy
private

ARKode solution vector.

Definition at line 878 of file arkode.h.

## ◆ abs_tolls

template<typename VectorType = Vector<double>>
 N_Vector SUNDIALS::ARKode< VectorType >::abs_tolls
private

ARKode absolute tolerances vector.

Definition at line 883 of file arkode.h.

## ◆ communicator

template<typename VectorType = Vector<double>>
 MPI_Comm SUNDIALS::ARKode< VectorType >::communicator
private

MPI communicator. SUNDIALS solver runs happily in parallel. Note that if the library is compiled without MPI support, MPI_Comm is aliased as int.

Definition at line 890 of file arkode.h.

## ◆ mem

template<typename VectorType = Vector<double>>
 private

Memory pool of vectors.

Definition at line 895 of file arkode.h.

The documentation for this class was generated from the following files: