 |
deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
|
Loading...
Searching...
No Matches
|
deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
|
#include <deal.II/sundials/arkode_stepper.h>
This class provides a wrapper to ARKStep time-stepping module. It is one of the possible steppers one can use together with the ARKode class.
ARKStep solves ODE initial value problems (IVPs) in \(R^N\). These problems should be posed in split, linearly-implicit 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:
ARKStep 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, ARKStep 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, ARKStep 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, ARKStep takes some shortcuts that allow a faster solution process.
For systems of either type, ARKStep 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 step 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, ARKStep 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, then jacobian_times_vector() should be supplied. If it 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 ARKStep 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 ARKStep, 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 ARKStep is provided in ARKStep documentation of ARKode.
The user has to provide the implementation of at least one (or both) of the following std::functions:
If the mass matrix is different from the identity, the user should supply
If the use of a Newton method is desired, then the user should also supply jacobian_times_vector(). jacobian_times_vector_setup() is optional.
A SUNDIALS default solver (SPGMR) is used to solve the linear systems. To use a custom linear solver for the mass matrix and/or Jacobian, set:
To use a custom preconditioner with either a default or custom linear solver, set:
Any other custom settings of the ARKStep object can be specified in
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{pmatrix} 0 & 1 \\ -k^2 &0 \end{pmatrix} \]
and \(y(0)=(0, k)^T\).
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:
Definition at line 401 of file arkode_stepper.h.
Classes | |
| class | AdditionalData |
Public Member Functions | |
| ARKStepper (const AdditionalData &data=AdditionalData()) | |
| void * | get_arkode_memory () const override |
Public Attributes | |
| std::function< void(const double t, const VectorType &y, VectorType &explicit_f)> | explicit_function |
| std::function< void(const double t, const VectorType &y, VectorType &res)> | implicit_function |
| std::function< void(const double t, const VectorType &v, VectorType &Mv)> | mass_times_vector |
| std::function< void(const double t)> | mass_times_vector_setup |
| std::function< void(const VectorType &v, VectorType &Jv, const double t, const VectorType &y, const VectorType &fy)> | jacobian_times_vector |
| std::function< void(const double t, const VectorType &y, const VectorType &fy)> | jacobian_times_vector_setup |
| LinearSolveFunction< VectorType > | solve_linearized_system |
| LinearSolveFunction< VectorType > | solve_mass |
| std::function< void(const double t, const VectorType &y, const VectorType &fy, const VectorType &r, VectorType &z, const double gamma, const double tol, const int lr)> | jacobian_preconditioner_solve |
| std::function< void(const double t, const VectorType &y, const VectorType &fy, const int jok, int &jcur, const double gamma)> | jacobian_preconditioner_setup |
| std::function< void(const double t, const VectorType &r, VectorType &z, const double tol, const int lr)> | mass_preconditioner_solve |
| std::function< void(const double t)> | mass_preconditioner_setup |
| std::function< void(void *arkode_mem)> | custom_setup |
Private Types | |
| using | CallbackContext = typename ARKodeStepper< VectorType >::template CallbackContext< ARKStepper< VectorType > > |
| using | ARKodeMemoryPtr = typename ARKodeStepper< VectorType >::ARKodeMemoryPtr |
Private Member Functions | |
| void | reinit (double t0, const VectorType &y0, internal::InvocationContext inv_ctx) override |
| void | setup_system_solver (const VectorType &solution, internal::InvocationContext inv_ctx) |
| void | setup_mass_solver (const VectorType &solution, internal::InvocationContext inv_ctx) |
Private Attributes | |
| AdditionalData | data |
| std::unique_ptr< internal::LinearSolverWrapper< VectorType > > | linear_solver |
| std::unique_ptr< internal::LinearSolverWrapper< VectorType > > | mass_solver |
| ARKodeMemoryPtr | arkode_mem |
| CallbackContext | callback_ctx |
|
private |
Definition at line 953 of file arkode_stepper.h.
|
private |
Definition at line 956 of file arkode_stepper.h.
| SUNDIALS::ARKStepper< VectorType >::ARKStepper | ( | const AdditionalData & | data = AdditionalData() | ) |
Constructor, with class parameters set by the AdditionalData object.
| data | ARKStep configuration data |
|
overridevirtual |
Provides user access to the internally used ARKODE memory.
This functionality is intended for users who wish to query additional information directly from the ARKODE integrator, refer to the ARKODE manual for the various ARKStepGet... functions. The ARKStepSet... functions should not be called since this might lead to conflicts with various settings that are performed by this ARKodeStepper object.
Implements SUNDIALS::ARKodeStepper< VectorType >.
|
overrideprivatevirtual |
Rebuild the stepper at a given time instance and for a given state vector. Required by the ARKodeStepper interface.
| t0 | Time instance that serves as starting time |
| y0 | Initial state vector whose layout is used for initialization of the internal ARKODE vectors |
| inv_ctx | Invocation context that provides access to the SUNContext object and the exception pointer managed by the caller |
Implements SUNDIALS::ARKodeStepper< VectorType >.
|
private |
Set up the (non)linear solver and preconditioners in the ARKODE memory object based on the user-specified functions.
| solution | The solution vector which is used as a template to create new vectors. |
| inv_ctx | Invocation context that provides access to the SUNContext object and the exception pointer managed by the caller |
|
private |
Set up the solver and preconditioner for a non-identity mass matrix in the ARKODE memory object based on the user-specified functions.
| solution | The solution vector which is used as a template to create new vectors. |
| inv_ctx | Invocation context that provides access to the SUNContext object and the exception pointer managed by the caller |
| std::function< void(const double t, const VectorType &y, VectorType &explicit_f)> SUNDIALS::ARKStepper< 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.
Definition at line 572 of file arkode_stepper.h.
| std::function<void(const double t, const VectorType &y, VectorType &res)> SUNDIALS::ARKStepper< 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.
Definition at line 591 of file arkode_stepper.h.
| std::function<void(const double t, const VectorType &v, VectorType &Mv)> SUNDIALS::ARKStepper< VectorType >::mass_times_vector |
A function object that users may supply and that is intended to compute the product of the mass matrix with a given vector v. This function will be called by ARKode (possibly several times) after mass_times_vector_setup() has been called at least once. ARKode tries to do its best to call mass_times_vector_setup() the minimum amount of times.
A call to this function should store in Mv the result of \(M\) applied to v.
Definition at line 612 of file arkode_stepper.h.
| std::function<void(const double t)> SUNDIALS::ARKStepper< VectorType >::mass_times_vector_setup |
A function object that users may supply and that is intended to set up 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 guaranteed to be called by ARKode at least once, before any call to mass_times_vector().
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 mass_times_vector() function, then the identity is used. If the mass_times_vector_setup() function is not provided, then mass_times_vector() 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 should be called to assemble the matrix. Subsequent calls (possibly more than one) to mass_times_vector() can assume that this function has been called at least once.
| t | The current evaluation time |
Definition at line 648 of file arkode_stepper.h.
| std::function<void(const VectorType &v, VectorType &Jv, const double t, const VectorType &y, const VectorType &fy)> SUNDIALS::ARKStepper< VectorType >::jacobian_times_vector |
A function object that users may supply and that is intended to compute the product of the Jacobian matrix with a given vector v. The Jacobian here refers to \(J=\frac{\partial f_I}{\partial y}\), i.e., the Jacobian of the user-specified implicit_function.
A call to this function should store in Jv the result of \(J\) applied to v.
Arguments to the function are
| [in] | v | The vector to be multiplied by the Jacobian |
| [out] | Jv | The vector to be filled with the product J*v |
| [in] | t | The current time |
| [in] | y | The current \(y\) vector for the current ARKode internal step |
| [in] | fy | The current value of the implicit right-hand side at y, \(f_I (t_n, y)\). |
Definition at line 681 of file arkode_stepper.h.
| std::function< void(const double t, const VectorType &y, const VectorType &fy)> SUNDIALS::ARKStepper< VectorType >::jacobian_times_vector_setup |
A function object that users may supply and that is intended to set up all data necessary for the application of jacobian_times_vector(). This function is called by ARKode any time a Jacobian update is required. The user should compute the Jacobian (or update all the variables that allow the application of Jacobian). This function is guaranteed to be called by ARKode at least once, before any call to jacobian_times_vector().
If the jacobian_times_vector_setup() function is not provided, then jacobian_times_vector() should do all the work by itself.
If the user uses a matrix based computation of the Jacobian, then this is the right place where an assembly routine should be called to assemble the matrix. Subsequent calls (possibly more than one) to jacobian_times_vector() can assume that this function has been called at least once.
| t | The current time |
| y | The current ARKode internal solution vector \(y\) |
| fy | The implicit right-hand side function evaluated at the current time \(t\) and state \(y\), i.e., \(f_I(y,t)\) |
Definition at line 720 of file arkode_stepper.h.
| LinearSolveFunction<VectorType> SUNDIALS::ARKStepper< VectorType >::solve_linearized_system |
A LinearSolveFunction object that users may supply and that is intended to solve the linearized system \(Ax=b\), where \(A = M-\gamma J\) is the Jacobian of the nonlinear residual. The application of the mass matrix \(M\) and Jacobian \(J\) are known through the functions mass_times_vector() and jacobian_times_vector() and \(\gamma\) is a factor provided by SUNDIALS. The matrix-vector product \(Ax\) is encoded in the supplied SundialsOperator. If a preconditioner was set through jacobian_preconditioner_solve(), it is encoded in the SundialsPreconditioner. If no preconditioner was supplied this way, the preconditioner is the identity matrix, i.e., no preconditioner. The user is free to use a custom preconditioner in this function object that is not supplied through SUNDIALS.
If you do not specify a solve_linearized_system() function, then a SUNDIALS packaged SPGMR solver with default settings is used.
For more details on the function type refer to LinearSolveFunction.
Definition at line 748 of file arkode_stepper.h.
| LinearSolveFunction<VectorType> SUNDIALS::ARKStepper< VectorType >::solve_mass |
A LinearSolveFunction object that users may supply and that is intended to solve the mass system \(Mx=b\). The matrix-vector product \(Mx\) is encoded in the supplied SundialsOperator. If a preconditioner was set through mass_preconditioner_solve(), it is encoded in the SundialsPreconditioner. If no preconditioner was supplied this way, the preconditioner is the identity matrix, i.e., no preconditioner. The user is free to use a custom preconditioner in this function object that is not supplied through SUNDIALS.
The user must specify this function if a non-identity mass matrix is used and applied in mass_times_vector().
For more details on the function type refer to LinearSolveFunction.
Definition at line 772 of file arkode_stepper.h.
| std::function<void(const double t, const VectorType &y, const VectorType &fy, const VectorType &r, VectorType &z, const double gamma, const double tol, const int lr)> SUNDIALS::ARKStepper< VectorType >::jacobian_preconditioner_solve |
A function object that users may supply to either pass a preconditioner to a SUNDIALS built-in solver or to apply a custom preconditioner within the user's own linear solve specified in solve_linearized_system().
This function should compute the solution to the preconditioner equation \(Pz=r\) and store it in z. In this equation \(P\) should approximate the Jacobian \(M-\gamma J\) of the nonlinear system.
| [in] | t | The current time |
| [in] | y | The current \(y\) vector for the current ARKode internal step |
| [in] | fy | The current value of the implicit right-hand side at y, \(f_I (t_n, y)\). |
| [in] | r | The right-hand side of the preconditioner equation |
| [out] | z | The solution of applying the preconditioner, i.e., solving \(Pz=r\) |
| [in] | gamma | The value \(\gamma\) in the preconditioner equation |
| [in] | tol | The tolerance up to which the system should be solved |
| [in] | lr | An input flag indicating whether the preconditioner solve is to use the left preconditioner (lr = 1) or the right preconditioner (lr = 2). Only relevant if used with a SUNDIALS packaged solver. If used with a custom solve_mass() function this will be set to zero. |
Definition at line 813 of file arkode_stepper.h.
| std::function<void(const double t, const VectorType &y, const VectorType &fy, const int jok, int &jcur, const double gamma)> SUNDIALS::ARKStepper< VectorType >::jacobian_preconditioner_setup |
A function object that users may supply to set up a preconditioner specified in jacobian_preconditioner_solve().
This function should prepare the solution of the preconditioner equation \(Pz=r\). In this equation \(P\) should approximate the Jacobian \(M-\gamma J\) of the nonlinear system.
If the jacobian_preconditioner_setup() function is not provided, then jacobian_preconditioner_solve() should do all the work by itself.
| [in] | t | The current time |
| [in] | y | The current \(y\) vector for the current ARKode internal step |
| [in] | fy | The current value of the implicit right-hand side at y, \(f_I (t_n, y)\). |
| [in] | jok | An input flag indicating whether the Jacobian-related data needs to be updated. The jok argument provides for the reuse of Jacobian data in the preconditioner solve function. When jok = SUNFALSE, the Jacobian-related data should be recomputed from scratch. When jok = SUNTRUE the Jacobian data, if saved from the previous call to this function, can be reused (with the current value of gamma). A call with jok = SUNTRUE can only occur after a call with jok = SUNFALSE. |
| [out] | jcur | On output this should be set to SUNTRUE if Jacobian data was recomputed, or set to SUNFALSE if Jacobian data was not recomputed, but saved data was still reused. |
| [in] | gamma | The value \(\gamma\) in \(M-\gamma J\). The preconditioner should approximate the inverse of this matrix. |
Definition at line 863 of file arkode_stepper.h.
| std::function<void(const double t, const VectorType &r, VectorType &z, const double tol, const int lr)> SUNDIALS::ARKStepper< VectorType >::mass_preconditioner_solve |
A function object that users may supply to either pass a preconditioner to a SUNDIALS built-in solver or to apply a custom preconditioner within the user's own linear solve specified in solve_mass().
This function should compute the solution to the preconditioner equation \(Pz=r\) and store it in z. In this equation \(P\) should approximate the mass matrix \(M\).
| [in] | t | The current time |
| [in] | r | The right-hand side of the preconditioner equation |
| [out] | z | The solution of applying the preconditioner, i.e., solving \(Pz=r\) |
| [in] | gamma | The value \(\gamma\) in the preconditioner equation |
| [in] | tol | The tolerance up to which the system should be solved |
| [in] | lr | An input flag indicating whether the preconditioner solve is to use the left preconditioner (lr = 1) or the right preconditioner (lr = 2). Only relevant if used with a SUNDIALS packaged solver. If used with a custom solve_mass() function this will be set to zero. |
Definition at line 897 of file arkode_stepper.h.
| std::function<void(const double t)> SUNDIALS::ARKStepper< VectorType >::mass_preconditioner_setup |
A function object that users may supply to set up a preconditioner specified in mass_preconditioner_solve().
This function should prepare the solution of the preconditioner equation \(Pz=r\). In this equation \(P\) should approximate the mass matrix \(M\).
If the mass_preconditioner_setup() function is not provided, then mass_preconditioner_solve() should do all the work by itself.
| [in] | t | The current time |
Definition at line 923 of file arkode_stepper.h.
| std::function<void(void *arkode_mem)> SUNDIALS::ARKStepper< VectorType >::custom_setup |
A function object that users may supply and which is intended to perform custom settings on the supplied arkode_mem object. Refer to the SUNDIALS documentation for valid options.
For instance, the following code specifies the fraction of the estimated explicitly stable step to use as 0.25 instead of the default value 0.5:
| arkode_mem | pointer to the ARKODE memory block which can be used for custom calls to ARKStepSet... methods. |
Definition at line 950 of file arkode_stepper.h.
|
private |
ARKStepper configuration data.
Definition at line 1003 of file arkode_stepper.h.
|
private |
Linear solver for applying the Jacobian of the implicit function in the nonlinear solver. This is used if the user provides both jacobian_times_vector() and solve_linearized_system().
Definition at line 1010 of file arkode_stepper.h.
|
private |
Linear solver for applying the mass matrix. This is used if the user provides both mass_times_vector() and solve_mass().
Definition at line 1016 of file arkode_stepper.h.
|
private |
ARKODE memory object. Declared after linear_solver and mass_solver so that it is destroyed first (C++ destroys members in reverse declaration order), ensuring ARKStepFree is called before the solver wrappers free their SUNLinearSolver objects.
Definition at line 1024 of file arkode_stepper.h.
|
private |
ARKStepper callback context.
Definition at line 1029 of file arkode_stepper.h.
1.9.8