|
KASKADE 7 development version
|
Spectral Deferred Correction methods for time stepping. More...
Classes | |
| class | Kaskade::SDCTimeGrid |
| Abstract base class of time grids for (block) spectral defect correction methods. More... | |
| class | Kaskade::SDCGridImplementationBase |
| A convenience class simplifying the implementation of different SDCTimeGrid derived classes. More... | |
| class | Kaskade::LobattoTimeGrid |
| spectral time grid for defect correction methods with Lobatto points on \( [a,b]\) More... | |
| class | Kaskade::RadauTimeGrid |
| spectral time grid for defect correction methods with Radau points on \( [a,b]\). More... | |
| class | Kaskade::GaussTimeGrid |
| spectral time grid for defect correction methods with Gauss points More... | |
Functions | |
| template<class Matrix , class Vectors , class ReactionDerivatives , class Solver > | |
| Matrix::field_type | Kaskade::sdcIterationStep (SDCTimeGrid const &grid, SDCTimeGrid::RealMatrix const &Shat, Solver const &solve, Matrix const &M, Matrix const &A, Vectors const &r, ReactionDerivatives const &f_u, Vectors const &u, Vectors &du) |
| A single spectral defect correction iteration sweep. More... | |
| template<class Matrix , class Vectors , class Solver > | |
| Matrix::field_type | Kaskade::sdcIMEXStep (SDCTimeGrid const &grid, SDCTimeGrid::RealMatrix const &Shat, Solver const &solve, Matrix const &M, Matrix const &A, Vectors const &r, Vectors const &u, Vectors &du) |
| A single IMEX SDIRK spectral deferred correction sweep. More... | |
| template<class Matrix , class Vectors , class ReactionDerivatives , class Solver > | |
| Matrix::field_type | Kaskade::sdcIterationStep2 (SDCTimeGrid const &grid, SDCTimeGrid::RealMatrix const &Shat, Solver const &solve, Matrix const &M, Matrix const &A, Vectors const &rUi, ReactionDerivatives const &f_u, Vectors const &Mdu, Vectors &du) |
| A single spectral defect correction iteration sweep. More... | |
| template<class Matrix , class Vectors , class ReactionDerivatives > | |
| double | Kaskade::waveformRelaxationStep2 (SDCTimeGrid const &grid, Matrix const &M, Matrix const &A, Vectors const &rUi, ReactionDerivatives const &f_u, Vectors &du) |
| A single waveform relaxation iteration. More... | |
Approximate integration matrices | |
The convergence properties of SDC are essentially defined by the selection of an approximate integration matrix \( \hat S \) acting as a lower triangular preconditioner for the spectral integration matrix \( S \). We provide several choices. | |
| void | Kaskade::eulerIntegrationMatrix (SDCTimeGrid const &grid, SDCTimeGrid::RealMatrix &Shat) |
| Triangular approximate integration matrix \( \hat S \) corresponding to the implicit Euler integrator. More... | |
| SDCTimeGrid::RealMatrix | Kaskade::eulerIntegrationMatrix (SDCTimeGrid const &grid) |
| Triangular approximate integration matrix \( \hat S \) corresponding to the implicit Euler integrator. More... | |
| void | Kaskade::luIntegrationMatrix (SDCTimeGrid const &grid, SDCTimeGrid::RealMatrix &Shat) |
| Triangular approximate integration matrix \( \hat S \) corresponding to a specialized DIRK sweep. More... | |
| SDCTimeGrid::RealMatrix | Kaskade::luIntegrationMatrix (SDCTimeGrid const &grid) |
| Triangular approximate integration matrix \( \hat S \) corresponding to a specialized DIRK sweep. More... | |
| SDCTimeGrid::RealMatrix | Kaskade::singleEulerIntegrationMatrix (SDCTimeGrid const &grid) |
| Triangular approximate integration matrix \( \hat S \) corresponding to an Euler sweep with constant "step size". More... | |
| SDCTimeGrid::RealMatrix | Kaskade::SDIRKIntegrationMatrix (SDCTimeGrid const &grid) |
| Triangular approximate integration matrix \( \hat S \) corresponding to an SDIRK sweep (i.e. with constant "step size"). More... | |
Spectral Deferred Correction methods for time stepping.
SDC is a class of flexible single step time integration methods that can be interpreted as
In particular the second aspect makes SDC methods attractive for combining them with inexact solves, mesh adaptivity, splitting methods, and multirate integration.
Kaskade 7 provides an implementation of SDC methods with a focus on parabolic problems. Several different interfaces are provided, tailored towards different problem settings and implementing different algorithmic variants.
| SDCTimeGrid::RealMatrix Kaskade::eulerIntegrationMatrix | ( | SDCTimeGrid const & | grid | ) |
Triangular approximate integration matrix \( \hat S \) corresponding to the implicit Euler integrator.
This is a convenience overload, see eulerIntegrationMatrix(SDCTimeGrid const& grid, SDCTimeGrid::RealMatrix& Shat) for full documentation.
| void Kaskade::eulerIntegrationMatrix | ( | SDCTimeGrid const & | grid, |
| SDCTimeGrid::RealMatrix & | Shat | ||
| ) |
Triangular approximate integration matrix \( \hat S \) corresponding to the implicit Euler integrator.
This computes a lower "triangular" matrix \( \hat{S} \in \mathbb{R}^{n\times n+1} \) with \( \hat{S}_{ij} = 0 \) for \( j > i+1 \) that can be used to formulate SDC sweeps.
This is \( \hat{S}_{i,i+1} = t_{i+1}-t_i \), all other entries zero. For quadrature, this is the right-looking rectangular rule. In an SDC sweep, this corresponds to an implicit Euler scheme for integrating the defect equation.
| [in] | grid | the collocation grid |
| [out] | Shat | the approximate integration matrix \( \hat S \) |
| SDCTimeGrid::RealMatrix Kaskade::luIntegrationMatrix | ( | SDCTimeGrid const & | grid | ) |
Triangular approximate integration matrix \( \hat S \) corresponding to a specialized DIRK sweep.
This is a convenience overload, see luIntegrationMatrix(SDCTimeGrid const& grid, SDCTimeGrid::RealMatrix& Shat) for full documentation.
| void Kaskade::luIntegrationMatrix | ( | SDCTimeGrid const & | grid, |
| SDCTimeGrid::RealMatrix & | Shat | ||
| ) |
Triangular approximate integration matrix \( \hat S \) corresponding to a specialized DIRK sweep.
This computes a lower "triangular" matrix \( \hat{S} \in \mathbb{R}^{n\times n+1} \) with \( \hat{S}_{ij} = 0 \) for \( j > i+1 \) that can be used to formulate SDC sweeps.
If \( S = [0 \, \tilde S] \) and \( \tilde S^T = LU \), this is \( \hat S = [0 \, U^T] \), which, in an SDC sweep corresponds to a cascade of tailor-made multi-step integrators forming one DIRK sweep. The consequence is that on the Dahlquist test equation \( \dot u = \lambda u \), in the limit \( \lambda \to \infty \) the SDC iteration matrix \( G = I - (U^T)^{-1} \tilde S \) is nilpotent and hence has a spectral radius of 0, which yields fast convergence for stiff components. Standard Euler-based SDC, in contrast, has a larger spectral radius for stiff components.
| [in] | grid | the collocation grid |
| [out] | Shat | the approximate integration matrix \( \hat S \) |
| Matrix::field_type Kaskade::sdcIMEXStep | ( | SDCTimeGrid const & | grid, |
| SDCTimeGrid::RealMatrix const & | Shat, | ||
| Solver const & | solve, | ||
| Matrix const & | M, | ||
| Matrix const & | A, | ||
| Vectors const & | r, | ||
| Vectors const & | u, | ||
| Vectors & | du | ||
| ) |
A single IMEX SDIRK spectral deferred correction sweep.
This is a specialization of the more general sdcIterationStep (see there for full documentation), but is more efficient (and convenient) for some classes of problems. In particular it
These two simplifications lead to the system matrices \( J_i = M-sA \) being identical, such that they neither need to be formed nor factorized (or preconditioner be constructed) over and over again. Consequently, the solver object is assumed to know the matrix \( J \) already.
If the time step remains constant over some steps, the same solver object (containing a factorization or a preconditioner) can be used unchanged for all those time steps.
| solve | a callable object that supports solve(x,b) for approximately solving \( J x = b \), where \( J=M-sA \) and \( x \) is the computed solution that on entry contains an initial guess that may or may not be used |
| Matrix::field_type Kaskade::sdcIterationStep | ( | SDCTimeGrid const & | grid, |
| SDCTimeGrid::RealMatrix const & | Shat, | ||
| Solver const & | solve, | ||
| Matrix const & | M, | ||
| Matrix const & | A, | ||
| Vectors const & | r, | ||
| ReactionDerivatives const & | f_u, | ||
| Vectors const & | u, | ||
| Vectors & | du | ||
| ) |
A single spectral defect correction iteration sweep.
SDC time stepping
This function performs one spectral defect correction (SDC) iteration for the abstract reaction diffusion equation
\[ M \dot u = A u + f(u) =: r(u) \]
using the linearly implicit Euler method on the given time grid. The solution interpolation vector u contains the approximate values of \( u \) at the time nodes \( t_i \) in a Lagrangian FE basis. Here, A and M are fixed matrices, such that locally we perform a method of lines in this time step.
SDC iterations are inexact Newton iterations for the Fredholm collocation time discretization of an ODE on a time grid \( t_0, \dots, t_n \):
\[ M (u_{i+1} - u_i) = \int_{\tau=t_i}^{t_{i+1}} p(\tau) \,d\tau \quad i=0,\dots,n-1, \]
where \( p(t_i) = r(u_i) = Au_i + f(u_i) \) represents the (probably polynomial) interpolant of the right hand side. By linearity, the integral can be expressed as a linear combination of the right hand side values:
\[ M (u_{i+1} - u_i) = \sum_{j=0}^n S_{ij} r(u_j) \]
Applying Newton's method to \( F(u) = M (u_{i+1} - u_i) - \sum_{j=0}^n S_{ij} r(u_j) \) yields \( F'(u^k)\delta u^k = -F(u^k), \quad \delta u_i^k = u_i^{k+1} - u_i^k \), or, more elaborate,
\[ M(\delta u_{i+1}^k - \delta u_i^k) - \sum_{j=0}^n S_{ij} r'(u_j^k) \delta u_j^k = -M (u_{i+1}^k - u_i^k) + \sum_{j=0}^n S_{ij} r(u_j^k) =: R_i, \quad i=0,\dots,n-1. \]
The time-global coupling makes this system difficult to solve. An approximation of the quadrature coefficients \( S_{ij} \) by a triangular \( \hat S_{ij} \) yields the simpler system
\[ M(\delta u_{i+1}^k - \delta u_i^k) - \hat{S}_{i,i+1} r'(u_{i+1}^k) \delta u_{i+1}^k = R_i + \sum_{j=0}^i \hat{S}_{i,j} r'(u_{j}^k) \delta u_j^k, \quad i=0,\dots,n-1, \]
starting with \( \delta u_0^k = 0 \), or, equivalently,
\[ (M - \hat{S}_{i,i+1} r'(u_{i+1}^k)) \delta u_{i+1}^k = M\delta u_i^k + R_i + \sum_{j=1}^i \hat{S}_{i,j} r'(u_{j}^k) \delta u_j^k, \quad i=0,\dots,n-1. \]
These \( n \) systems with matrix \( J_i = M - \hat{S}_{i,i+1} r'(u_{i+1}^k) \) need to be solved, for which a solver needs to be provided.
The norm of the correction \( \delta u^k \) is returned, i.e.
\[ \left( (t_n-t_0)^{-1}\sum_{i=0}^{n-1} (t_{i+1}-t_i)(\delta u_i^k)^T (M- \hat{S}_{i,i+1} r'(u_{i+1}^k)) \delta u_i^k \right)^{1/2}. \]
Reaction derivatives
If the nonlinear reaction \( f \) is a Nemyckii operator and thus acts locally (which is quite common) and the FE discretization is done by Lagrange finite elements (also quite common), the derivative \( f' \) can be well approximated by a diagonal matrix, which simplifies the derivative computation and makes the construction of \( J_i \) faster.
| Matrix | a sparse matrix type, usually Dune::BCRSMatrix or a NumaBCRSMatrix |
| Vectors | a container type with elements from the domain type of Matrix |
| Solver | a callable type for (approximately) solving linear equation systems |
| ReactionDerivatives | a container type with sparse matrix elements (usually BCRSMatrix or NumaBCRSMatrix) |
| [in] | grid | the collocation time grid with \( n+1 \) points |
| [in] | Shat | the triangular quadrature matrix approximation \( \hat S \in \mathbb{R}^{n\times n+1} \) (see eulerIntegrationMatrix()) |
| [in] | solve | a callable that supports solve(J,x,b) giving an approximative solution of \( Jx = b \), where J is of type Matrix and \( x \) is the computed solution that on entry contains an initial guess that may or may not be used |
| [in] | M | the mass matrix |
| [in] | A | the stiffness matrix \( A \) (usually negative semidefinite) |
| [in] | r | the right hand sides \( r_i = Au_i +f(u_i)\) at the collocation times (including interval start and end point), size n+1 |
| [in] | f_u | the reaction term derivatives (a collection of sparse matrices, the sparsity pattern of which is a subset of that of A and M), size n+1 |
| [in] | u | the current iterate (values of \( u \) at the collocation times including start and end point), size n+1 |
| [out] | du | the approximate Newton correction \( \delta u \), size >=n+1 |
| Matrix::field_type Kaskade::sdcIterationStep2 | ( | SDCTimeGrid const & | grid, |
| SDCTimeGrid::RealMatrix const & | Shat, | ||
| Solver const & | solve, | ||
| Matrix const & | M, | ||
| Matrix const & | A, | ||
| Vectors const & | rUi, | ||
| ReactionDerivatives const & | f_u, | ||
| Vectors const & | Mdu, | ||
| Vectors & | du | ||
| ) |
A single spectral defect correction iteration sweep.
This is an overload with slightly different interface (required for some more or less arcane algorithmic variants). The only interface difference is that instead of values \( u_i \), the actually required differences \( M (u_{i+1}-u_i) \) are provided. Those are computed explicitly in sdcIterationStep.
| Matrix | a sparse matrix type, usually Dune::BCRSMatrix or a NumaBCRSMatrix |
| Vectors | a container type with elements from the domain type of Matrix |
| ReactionDerivatives | either (i) a container type with sparse matrix elements (usually BCRSMatrix or NumaBCRSMatrix) or (ii) nullptr_t, signalling that reaction shall be treated explicitly (i.e. as if \( f' = 0 \). |
| Solver |
| [in] | grid | the collocation time grid with n+1 points |
| [in] | Shat | the triangular quadrature matrix approximation \( \hat S \) |
| [in] | solve | a callable that supports solve(A,x,b) giving an approximative solution of \( Ax = b \), where A is of type Matrix |
| [in] | M | the mass matrix |
| [in] | A | the stiffness matrix (usually negative semidefinite) |
| [in] | rUi | the right hand sides at the collocation times (including interval start and end point), size n+1 |
| [in] | f_u | the reaction term derivatives (a collection of sparse matrices, the sparsity pattern of which is a subset of that of A and M), size n+1 |
| [in] | Mdu | the current iterate differences (values of \( M(u_{i}-u_{i+1}) \) for \( i=0,\dots,n-1 \)) |
| [out] | du | the approximate Newton correction \( \delta u \) |
Definition at line 716 of file sdc.hh.
Referenced by Kaskade::sdcIMEXStep(), and Kaskade::sdcIterationStep().
| SDCTimeGrid::RealMatrix Kaskade::SDIRKIntegrationMatrix | ( | SDCTimeGrid const & | grid | ) |
Triangular approximate integration matrix \( \hat S \) corresponding to an SDIRK sweep (i.e. with constant "step size").
The integration matrix has the same value \( s \) on all diagonal entries, such that the "Euler" steps in the sweep all read \( (M-sA) \delta u_i = \tau_i (\dots) \). This allows re-using a factorization of (or a preconditioner for) \( M-sA\) for all steps in the sweep.
On Radau grids, the worst case contraction rate for Dahlquist's test equation is slightly typically much better than that of the standard Euler sweep. This is important for problems where either mid- to high-frequent error terms are constantly excited (such as in nonlinear boundary conditions, e.g.), or where mid- to high-frequent errors spill over into error components that are conserved (such as the position error in monodomain equations). In problems like the linear heat equation, the difference to the single Euler sweep is barely noticeable, since the error components that the Euler sweep does not extinguish are eliminated in subsequent time steps.
For some background information, see Computation of integration matrices for SDC.
| SDCTimeGrid::RealMatrix Kaskade::singleEulerIntegrationMatrix | ( | SDCTimeGrid const & | grid | ) |
Triangular approximate integration matrix \( \hat S \) corresponding to an Euler sweep with constant "step size".
The integration matrix has the same value \( s \) on all diagonal entries, such that the Euler steps in the sweep all read \( (M-sA) \delta u_i = \tau_i (\dots) \). This allows re-using a factorization of (or a preconditioner for) \( M-sA\) for all steps in the sweep.
On Radau grids, the worst case contraction rate for Dahlquist's test equation is slightly worse than the pure Euler sweep for less than four collocation points, and slightly better for more collocation points.
For some background information, see Computation of integration matrices for SDC.
| double Kaskade::waveformRelaxationStep2 | ( | SDCTimeGrid const & | grid, |
| Matrix const & | M, | ||
| Matrix const & | A, | ||
| Vectors const & | rUi, | ||
| ReactionDerivatives const & | f_u, | ||
| Vectors & | du | ||
| ) |
A single waveform relaxation iteration.
This function performs one linearized waveform relaxation iteration for
\[ M \dot u = A u + M f(u) \]
using the Jacobi method. It approximately solves the linear defect system
\[ M \dot{\delta u} = A \delta u + M f'(u)\delta u + r, \quad r=Au+f(u)-\dot u \]
for the correction \( \delta u \) to be added to the provided approximation \( u \). The time evolution is discretized by a collocation method on the given time grid. Thus, for each entry \( u_i \), the collocation solution of the scalar ODE
\[ M_{ii}\dot{\delta u_i} = A_{ii} \delta u_i + M_{ii}f'(u_i) \delta u_i + r_i. \]
The solution interpolation vector
Linearized waveform relaxation works here as
The norm of the correction \( \delta u^k \) is returned, i.e.
\[ \]
| Matrix | a sparse matrix type, usually Dune::BCRSMatrix or a NumaBCRSMatrix |
| Vectors | a container type with elements from the domain type of Matrix |
| ReactionDerivatives |
| [in] | grid | the collocation time grid |
| [in] | Shat | the triangular quadrature matrix approximation \( \hat S \) |
| [in] | solve | a callable that supports solve(A,x,b) giving an approximative solution of \( Ax = b \), where A is of type Matrix |
| [in] | M | the mass matrix |
| [in] | A | the stiffness matrix (usually negative semidefinite) |
| [in] | rUi | the right hand sides at the collocation times |
| [in] | f_u | the reaction term derivatives (a container of sparse matrices, the sparsity pattern of which is a subset of that of A and M), size n+1 |
| [in] | Mdu | the current iterate differences (values of \( M(u_{i}-u_{i+1}) \) for \( i=0,\dots,n-1 \)) |
| [out] | du | the approximate Newton correction \( \delta u \) |
1.9.4