## 1 Introduction

Explicit time integrators for parabolic PDE are subject to the restrictive timestep limit , so A-stable integrators are essential. It is well known that although there are no A-stable explicit linear multistep methods and implicit multistep methods cannot be A-stable beyond order two, there exist A-stable and L-stable implicit Runge–Kutta (IRK) methods at all orders [hairer1993solving, wanner1996solving]. IRK methods offer an appealing combination of stability and high order; however, these methods are not widely used for PDE because they lead to large, strongly coupled linear systems. An -stage IRK system has

times as many degrees of freedom as the systems resulting from backward Euler or implicit trapezoidal rule discretization applied to the same equation set. Order-optimal preconditioners for such systems have been investigated in a series of papers

[mardal2007order, staff2006preconditioning]. In this paper we introduce a new block preconditioner for IRK methods, based on an LDU factorization. Solves on individual blocks are accomplished using a multigrid algorithm. We demonstrate the effectiveness of this preconditioner on the heat equation as a simple test problem, and find that it is scalable (-independent) and yields better timing results than other preconditioners currently in the literature.## 2 Implicit Runge–Kutta Methods

IRK methods are described in detail in, for example, [wanner1996solving]. The general -stage IRK method for requires the solution of simultaneous equations

for the stage variables , . The stage variables are then used to update the solution ,

The coefficients that define a given method are summarized in a Butcher table:

We list some common IRK methods with their order and stability properties in Table 1. As our goal is to use higher order L-stable methods, throughout this paper we consider Radau IIA and Lobatto IIIC methods. Of the common L-stable methods, these provide the highest order for a given number of stages, that is, for a given cost.

IRK Methods | Order | Stability |
---|---|---|

Gauss–Legendre (s) | A-stable | |

Radau IA (s) | L-stable | |

Radau IIA (s) | L-stable | |

Lobatto IIIA (s) | A-stable | |

Lobatto IIIC (s) | L-stable | |

Miller DIRK (2) | 2 | L-stable |

Miller DIRK (3) | 2 | L-stable |

Crouzeix SDIRK(2) | 2 | L-stable |

Crouzeix SDIRK(3) | 3 | L-stable |

## 3 Formulation of Test Problem

As a simple test problem, we consider the heat equation on some spatial domain over the time interval with Dirichlet boundary conditions:

(1) |

We first discretize in time using an IRK method. We write for and introduce stage variables for to . In strong form, the stage equations for an -stage IRK method are:

We discretize this system of time-independent PDE using the finite element method. Adopting test functions , we obtain a weak form as usual by multiplying the -th stage equation by the test function , integrating, and applying Green’s identity. This leads to the equations:

which must hold for all .
Note that the boundary terms vanish due to the boundary restriction of
to .
Therefore we need to find such that

Converting to weak form and discretizing with finite element basis functions , we get the following linear system, which must be solved for the stage variables :

(2) |

where the matrices and are given by

Thus we have an linear system to solve at each time step; this is the system that we need to precondition.

## 4 Preconditioning IRK methods applied to parabolic PDEs

In this paper, we develop a block upper triangular preconditioner and a block lower triangular preconditioner for the system (2).

Let be the -stage IRK coefficient matrix, that is, the matrix comprising the elements in the upper-right block of the Butcher table for the given IRK method.

Given an matrix and a matrix , recall that their Kronecker product is given by

Using this notation, we can write the block matrix (2) arising from the FEM discretization as

where is the identity matrix.

In this paper, we consider several block preconditioners of the general form

(3) |

where is either a diagonal or triangular matrix referred as the preconditioner coefficient matrix. We can, for example, define a block Jacobi preconditioner and a block Gauss–Seidel preconditioner using this notation as follows:

(4) | ||||

where is the diagonal part of and is the lower triangular part of . These preconditioners and their order optimality have been extensively studied in Mardal et al. [mardal2007order] and in G.A. Staff et al. [staff2006preconditioning]. We now introduce two new -based block triangular preconditioners for .

### 4.1 Ldu-based block triangular preconditioners

It has been shown in [staff2006preconditioning] that if all blocks of are well preconditioned by a preconditioner of the form (3), then all of the eigenvalues of the preconditioned system will be clustered, and the condition number of the preconditioned system can be approximated by , where is the preconditioner coefficient matrix. In other words, if we make an effective preconditioner for the Butcher coefficient matrix , then can be expected to be an effective (left) preconditioner for the system . (In practice, we use the same motivation in developing left and right preconditioners.) We begin, therefore, by preconditioning .

For the development of our preconditioners, we assume that the IRK coefficient matrix is invertible and has nonsingular leading principal submatrices. It is proven in [hairer1993solving] that the Butcher coefficient matrices for Gauss–Legendre, Radau IA, Radau IIA, and Lobatto IIIC are nonsingular for all . We have confirmed that the leading principal submatrices are all nonsingular for the coefficient matrices arising from Gauss–Legendre, Lobatto IIIC, Radau IA, and Radau IIA IRK methods for all stages that we use in this paper: through . Given these conditions, can be factored into without pivoting. We compute formally the factorization without pivoting

and consider and as preconditioners for . Since , we have

and | |||

The eigenvalues of and are all one, thus we expect to be a good right preconditioner and to be a good left preconditioner for . Since the eigenvalues (though not the condition numbers) of the product of two matrices do not depend on the order of the product, we consider both matrices as left and right preconditioners in our numerical experiments.

Using the condition number motivation from [staff2006preconditioning], we now define block upper and lower triangular preconditioners for of the form:

(5) | ||||

### 4.2 Application of the Preconditioners

Applying the block Jacobi preconditioner

to a vector

is straightforward, and for an -stage IRK scheme involves subsolves. The subsolve for the th block has the formIn practice and for efficiency, we use a single algebraic multigrid (AMG) V-cycle for each subsolve for the block .

Application of the block upper triangular preconditioner is done via back substitution and again involves subsolves. The subsolve on the th block in this case is of the form

In practice, we again use a single AMG V-cycle for each subsolve.

Application of the block lower triangular preconditioners and is done via forward substitution and again involves subsolves. The subsolves on the th block in this case are of the form

As with the previous preconditioners, we use a single AMG V-cycle for each subsolve.

### 4.3 Analysis of the Preconditioners

The preconditioners have been chosen because of their intended effect on the condition number of the preconditioned system. In Table 2, we examine the 2-norm condition number for the unpreconditioned system , and for preconditioned on the right by , , , and . The system is from a 2D problem using an -stage Radau IIA IRK method for ranging from 2 to 7. In all of these results, we have chosen and , where is the degree of the Lagrange polynomial basis functions in space. The preconditioners are constructed exactly. All of the preconditioners significantly reduce the condition number, with and giving the most improvement. gives the greatest improvement on the condition number for the lower-order stages and , and gives the greatest improvement for the higher-order stages through .

s | |||||
---|---|---|---|---|---|

2 | 240.37 | 3.23 | 1.75 | 5.32 | 2.48 |

3 | 502.53 | 5.66 | 2.58 | 11.18 | 2.66 |

4 | 746.23 | 8.54 | 3.63 | 18.23 | 3.04 |

5 | 959.16 | 11.76 | 5.08 | 26.53 | 3.21 |

6 | 1137.24 | 15.23 | 7.13 | 35.97 | 3.50 |

7 | 1281.47 | 18.90 | 10.05 | 46.48 | 3.67 |

With GMRES, the condition number does not always tell the whole story, and it is also important to consider how well the preconditioner clusters the eigenvalues. In Figure 1 we show eigenvalues plotted in the complex plane for the matrix and the same system right preconditioned with the same four preconditioners: , , , and . Here we use the same 2D problem as before, using an -stage Radau IIA IRK method for ranging from 2 to 5. We have again chosen and , where is the degree of the Lagrange polynomial basis functions in space. The preconditioners are constructed exactly. The original matrix has a number of eigenvalues very near zero. All of the preconditioners succeed in clustering the eigenvalues farther away from zero. For , and give the tightest clustering and have very similar eigenvalues. But as increases, maintains tighter clustering than . This fact is reflected in our numerical results presented in the next section.

### 4.4 Comparison between with optimal coefficients and preconditioners

In G.A. Staff et al. [staff2006preconditioning], the authors develop an optimal block lower triangular Gauss–Seidel preconditioner. This block Gauss–Seidel preconditioner is derived such that the preconditioner coefficients are computed from the following optimization problem

(6) |

That is, their optimal preconditioner, which we will refer to as , uses coefficients that optimize the condition number of the (left) preconditioned system.

In Table 3, we examine the 2-norm condition number for left preconditioned by and for 1D and 2D problems. We fix and , where is the degree of the Lagrange polynomial basis functions in space and is the order of the corresponding IRK method. We use Radau IIA for through and Lobatto IIIC for through . For the Radau IIA problems, the preconditioner yields a lower condition number than the optimized preconditioner for all stages . For the Lobatto IIIC problems, we have the opposite result with consistently having lower condition number than for all stages .

1D | 2D | |||

R IIA 2 | 1.58 | 1.26 | 1.60 | 1.26 |

R IIA 3 | 1.94 | 1.50 | 1.93 | 1.52 |

R IIA 4 | 2.18 | 1.74 | 2.17 | 1.77 |

R IIA 5 | 2.37 | 1.95 | 2.40 | 1.98 |

R IIA 6 | 3.01 | 2.14 | 3.18 | 2.18 |

L IIIC 2 | 1.48 | 2.68 | 1.54 | 2.78 |

L IIIC 3 | 3.20 | 6.68 | 3.32 | 6.94 |

L IIIC 4 | 4.92 | 10.43 | 5.24 | 10.93 |

Again we note that for GMRES, a lower condition number does not always indicate a superior preconditioner. In Figures 2 and 3 we see that more effectively clusters the eigenvalues for both Radau IIA and Lobatto IIIC for stages and .

## 5 Numerical Experiments

To test the performance of our preconditioners , , , and , we consider a problem (1) with domain and discretize using bivariate quadratic finite elements in space and -stage RadauIIA and Lobatto IIIC methods in time. Spatial discretization is done using the Sundance finite element toolkit [long_unified_2010]. The method of manufactured solutions is used to produce an exact solution to each linear subproblem so that we can check the error in addition to the residual.

A common practice in analysis of methods for spatiotemporal problems is to choose and increase the polynomial order of spatial discretization in lockstep with the order of time discretization. Here, we use a different approach. We use quadratic elements () throughout, and then for each we choose a timestep such that the spatial discretization error and temporal global truncation error are comparable. With a Radau -stage method we set so that the timestep is .

By holding the spatial order fixed, we are simulating conditions in which a modeler uses a high order integrator to obtain results of a specified accuracy with larger timestep and smaller computational cost.

We solve the linear system (2) using preconditioned GMRES. We construct our preconditioners as described in Section 4 using one AMG V-cycle for each subsolve. We use AMG from the IFISS software package [ifiss]. All results are computed in MATLAB on a machine with an Intel Core i7 1.80 GHz (Turbo Boost up to 4.00 GHz) processor and 8.00 GB RAM.

### 5.1 Comparison results between , , , and preconditioners

In Table 4 we report iteration counts and timing results for right-preconditioned GMRES to converge with relative tolerance for varying mesh sizes and varying number of stages (varying order). Throughout, all of the triangular preconditioners outperform the block Jacobi preconditioner , as expected. The -based block upper triangular preconditioner and the block Gauss–Seidel preconditioner perform similarly, with performing slightly better, especially for the higher stage methods. However, overall, the -based block lower triangular preconditioner performs better than all of the other preconditioners in both iteration count and timing. For 2-stage Radau IIA, and have similar performance, with performing slightly better on the smaller problems and performing slightly better on the larger problems. But as the number of stages increases (that is, with increasing order), the performance of improves over the other preconditioners giving consistently lower iteration counts and timing. For a given number of stages, all of the triangular preconditioners scale well with problem size, exhibiting little to no growth in iteration count with increasing problem size. However the preconditioner has the least increase in cost as the number of stages increases.

stage | DOF | |||||
---|---|---|---|---|---|---|

s=2 | 8 | 450 | 14 (0.03) | 8 (0.01) | 7 (0.01) | 7 (0.01) |

16 | 1922 | 14 (0.05) | 8 (0.02) | 7 (0.02) | 7 (0.02) | |

32 | 7938 | 14 (0.10) | 8 (0.06) | 7 (0.06) | 7 (0.05) | |

64 | 32,258 | 14 (0.40) | 7 (0.19) | 7 (0.19) | 7 (0.19) | |

128 | 130,050 | 14 (1.48) | 7 (0.80) | 7 (0.83) | 7 (0.82) | |

s=3 | 8 | 675 | 23 (0.06) | 10 (0.05) | 10 (0.03) | 9 (0.03) |

16 | 2883 | 22 (0.06) | 10 (0.03) | 10 (0.03) | 8 (0.02) | |

32 | 11,907 | 21 (0.17) | 10 (0.09) | 10 (0.08) | 8 (0.07) | |

64 | 48,387 | 22 (0.67) | 9 (0.29) | 10 (0.32) | 8 (0.26) | |

128 | 195,075 | 21 (2.90) | 8 (1.12) | 9 (1.25) | 8 (1.12) | |

s=4 | 8 | 900 | 32 (0.10) | 13 (0.08) | 13 (0.07) | 10 (0.04) |

16 | 3844 | 31 (0.11) | 12 (0.05) | 13 (0.05) | 10 (0.04) | |

32 | 15,876 | 30 (0.31) | 12 (0.13) | 12 (0.13) | 10 (0.12) | |

64 | 64,516 | 29 (1.16) | 11 (0.47) | 12 (0.52) | 9 (0.39) | |

128 | 260,100 | 27 (5.24) | 11 (2.10) | 12 (2.29) | 9 (1.73) | |

s=5 | 8 | 1125 | 42 (0.10) | 15 (0.08) | 16 (0.08) | 11 (0.06) |

16 | 4805 | 40 (0.23) | 15 (0.08) | 16 (0.08) | 11 (0.06) | |

32 | 19,845 | 39 (0.50) | 13 (0.18) | 15 (0.20) | 11 (0.17) | |

64 | 80,645 | 35 (1.80) | 13 (0.70) | 15 (0.81) | 11 (0.61) | |

128 | 325,125 | 33 (8.44) | 12 (2.98) | 14 (3.49) | 11 (2.75) | |

s=6 | 8 | 1350 | 53 (0.15) | 18 (0.08) | 19 (0.08) | 12 (0.04) |

16 | 5766 | 50 (0.31) | 17 (0.11) | 18 (0.11) | 12 (0.08) | |

32 | 23,814 | 46 (0.73) | 16 (0.27) | 18 (0.29) | 12 (0.19) | |

64 | 96,774 | 42 (2.66) | 15 (1.02) | 17 (1.14) | 12 (0.81) | |

128 | 390,150 | 38 (12.13) | 14 (4.27) | 17 (5.41) | 12 (3.56) | |

s=7 | 8 | 1575 | 62 (0.22) | 21 (0.10) | 23 (0.10) | 13 (0.09) |

16 | 6727 | 57 (0.40) | 20 (0.15) | 22 (0.16) | 13 (0.09) | |

32 | 27,783 | 52 (0.98) | 19 (0.37) | 21 (0.40) | 13 (0.26) | |

64 | 112,903 | 48 (3.64) | 18 (1.49) | 20 (1.61) | 12 (0.98) | |

128 | 455,175 | 44 (17.06) | 17 (6.29) | 20 (7.24) | 12 (4.42) |

We have also tested all four preconditioners for Lobatto IIIC methods. Table 5 shows iteration counts and timing results for right-preconditioned GMRES to converge with a relative tolerance of for the problem (1) using Lobatto IIIC methods in time and quadratic finite elements in space. For the Lobatto IIIC methods, the preconditioner outperforms all of the other preconditioners both in iteration count and timing, with more significant improvements for larger problems and larger number of stages (higher order).

stage | DOF | |||||
---|---|---|---|---|---|---|

s=2 | 8 | 450 | 17 (0.01) | 10 (0.01) | 7 (0.03) | 7 (0.01) |

16 | 1922 | 18 (0.03) | 10 (0.02) | 7 (0.02) | 8 (0.02) | |

32 | 7938 | 19 (0.13) | 10 (0.07) | 8 (0.06) | 8 (0.06) | |

64 | 32,258 | 19 (0.47) | 10 (0.26) | 7 (0.19) | 8 (0.22) | |

128 | 130,050 | 18 (2.19) | 10 (1.21) | 7 (0.86) | 8 (1.00) | |

s=3 | 8 | 675 | 27 (0.05) | 12 (0.03) | 11 (0.02) | 10 (0.02) |

16 | 2883 | 29 (0.07) | 12 (0.03) | 11 (0.03) | 10 (0.03) | |

32 | 11,907 | 27 (0.22) | 12 (0.10) | 11 (0.10) | 10 (0.08) | |

64 | 48,387 | 29 (0.96) | 11 (0.36) | 11 (0.34) | 10 (0.35) | |

128 | 195,075 | 27 (3.89) | 11 (1.54) | 10 (1.41) | 9 (1.27) | |

s=4 | 8 | 900 | 40 (0.07) | 16 (0.05) | 15 (0.05) | 12 (0.05) |

16 | 3844 | 40 (0.14) | 15 (0.06) | 15 (0.06) | 12 (0.05) | |

32 | 15,876 | 39 (0.41) | 14 (0.15) | 14 (0.15) | 12 (0.13) | |

64 | 64,516 | 38 (1.58) | 14 (0.60) | 13 (0.57) | 11 (0.47) | |

128 | 260,100 | 37 (7.58) | 13 (2.46) | 13 (2.44) | 11 (2.06) | |

s=5 | 8 | 1125 | 55 (0.08) | 19 (0.03) | 19 (0.03) | 13 (0.02) |

16 | 4805 | 55 (0.26) | 18 (0.09) | 18 (0.08) | 13 (0.07) | |

32 | 19,845 | 51 (0.57) | 17 (0.20) | 17 (0.22) | 13 (0.16) | |

64 | 80,645 | 49 (2.61) | 16 (0.85) | 16 (0.84) | 12 (0.64) | |

128 | 325,125 | 45 (12.01) | 15 (3.64) | 15 (3.69) | 12 (2.96) |

### 5.2 Comparison results between with optimal coefficients and preconditioners

In this section we compare our -based block lower triangular preconditioner with the optimal block lower triangular Gauss–Seidel preconditioner introduced in [staff2006preconditioning]. was developed using coefficients optimized to reduce the condition number of the left-preconditioned system .

In Table 6, we present results similar to those in Table 4. The problem set up is the same, but here we use left preconditioning for all of the preconditioners since the coefficients of were optimized for left preconditioning, and we include an additional column for . The first things we note is that although the idea of optimizing coefficients to reduce the condition number of the preconditioned system is a perfectly sensible idea, for this problem the unoptimized performs better than the optimized except on the smallest problems. For all numbers of stages, performs better than for problems sizes and . We observe in this table that our has the best performance overall. Although it performs slightly worse as a left preconditioner than it did in Table 4 as a right preconditioner, it nonetheless achieves the lowest iteration count and lowest timing of all of the preconditioners for all numbers of stages and all problem sizes.

stage | DOF | ||||||
---|---|---|---|---|---|---|---|

s=2 | 8 | 450 | 14 (0.02) | 9 (0.02) | 9 (0.01) | 7 (0.01) | 7 (0.01) |

16 | 1922 | 16 (0.03) | 9 (0.02) | 10 (0.02) | 8 (0.02) | 7 (0.02) | |

32 | 7938 | 16 (0.10) | 9 (0.06) | 10 (0.06) | 8 (0.05) | 7 (0.05) | |

64 | 32,258 | 17 (0.38) | 9 (0.24) | 10 (0.27) | 8 (0.23) | 7 (0.21) | |

128 | 130,050 | 18 (2.01) | 9 (1.16) | 10 (1.24) | 8 (1.03) | 7 (0.94) | |

s=3 | 8 | 675 | 23 (0.06) | 11 (0.04) | 13 (0.03) | 10 (0.04) | 9 (0.03) |

16 | 2883 | 25 (0.05) | 12 (0.03) | 13 (0.03) | 11 (0.03) | 9 (0.02) | |

32 | 11,907 | 25 (0.20) | 12 (0.09) | 14 (0.10) | 11 (0.08) | 9 (0.07) | |

64 | 48,387 | 27 (0.77) | 12 (0.36) | 14 (0.43) | 11 (0.34) | 9 (0.31) | |

128 | 195,075 | 28 (3.65) | 12 (1.65) | 14 (1.91) | 11 (1.56) | 9 (1.32) | |

s=4 | 8 | 900 | 32 (0.07) | 15 (0.06) | 18 (0.05) | 14 (0.06) | 11 (0.04) |

16 | 3844 | 34 (0.13) | 15 (0.06) | 19 (0.07) | 14 (0.06) | 11 (0.05) | |

32 | 15,876 | 36 (0.38) | 15 (0.17) | 19 (0.21) | 14 (0.16) | 11 (0.13) | |

64 | 64,516 | 36 (1.52) | 15 (0.70) | 19 (0.84) | 14 (0.64) | 11 (0.51) | |

128 | 260,100 | 36 (6.51) | 15 (2.88) | 19 (3.63) | 14 (2.74) | 11 (2.20) | |

s=5 | 8 | 1125 | 44 (0.11) | 18 (0.08) | 22 (0.06) | 17 (0.07) | 12 (0.06) |

16 | 4805 | 45 (0.21) | 18 (0.08) | 22 (0.10) | 17 (0.08) | 12 (0.06) | |

32 | 19,845 | 46 (0.61) | 18 (0.26) | 22 (0.28) | 17 (0.23) | 12 (0.17) | |

64 | 80,645 | 47 (2.21) | 18 (0.90) | 22 (1.18) | 17 (0.93) | 12 (0.74) | |

128 | 325,125 | 47 (11.15) | 18 (4.41) | 22 (5.45) | 17 (4.26) | 12 (3.00) | |

s=6 | 8 | 1350 | 55 (0.14) | 22 (0.08) | 24 (0.07) | 20 (0.07) | 13 (0.04) |

16 | 5766 | 55 (0.32) | 22 (0.14) | 24 (0.16) | 21 (0.14) | 14 (0.09) | |

32 | 23,814 | 56 (0.93) | 22 (0.38) | 24 (0.40) | 20 (0.34) | 14 (0.24) | |

64 | 96,774 | 56 (3.77) | 21 (1.49) | 24 (1.69) | 20 (1.42) | 13 (0.94) | |

128 | 390,150 | 56 (15.89) | 21 (6.21) | 23 (6.89) | 20 (5.63) | 13 (3.98) |

In Table 7, we present results similar to those in Table 6 but using Lobatto IIIC methods. The problem set up is again the same, using left preconditioning for all of the preconditioners. The results are very similar to those in the previous table. The unoptimized performs better than the optimized on all of these problems. Our has the best performance overall.

stage | DOF | ||||||
---|---|---|---|---|---|---|---|

s=2 | 8 | 450 | 18 (0.02) | 11 (0.02) | 18 (0.02) | 7 (0.02) | 8 (0.02) |

16 | 1922 | 20 (0.04) | 12 (0.02) | 20 (0.03) | 8 (0.02) | 8 (0.02) | |

32 | 7938 | 21 (0.15) | 12 (0.09) | 21 (0.14) | 8 (0.06) | 8 (0.06) | |

64 | 32,258 | 21 (0.56) | 12 (0.34) | 21 (0.56) | 8 (0.24) | 8 (0.23) | |

128 | 130,050 | 21 (2.40) | 12 (1.49) | 21 (2.50) | 8 (1.05) | 8 (1.05) | |

s=3 | 8 | 675 | 27 (0.06) | 14 (0.04) | 25 (0.04) | 11 (0.03) | 11 (0.03) |

16 | 2883 | 32 (0.08) | 14 (0.04) | 27 (0.06) | 12 (0.03) | 11 (0.03) | |

32 | 11,907 | 33 (0.23) | 15 (0.11) | 28 (0.20) | 12 (0.09) | 11 (0.09) | |

64 | 48,387 | 35 (1.06) | 15 (0.49) | 30 (0.98) | 12 (0.40) | 11 (0.37) | |

128 | 195,075 | 37 (3.89) | 15 (1.87) | 30 (3.76) | 12 (1.67) | 11 (1.58) | |

s=4 | 8 | 900 | 41 (0.09) | 18 (0.06) | 31 (0.07) | 15 (0.06) | 13 (0.05) |

16 | 3844 | 44 (0.13) | 18 (0.06) | 35 (0.11) | 16 (0.06) | 13 (0.05) | |

32 | 15,876 | 46 (0.51) | 18 (0.20) | 35 (0.38) | 16 (0.18) | 13 (0.15) | |

64 | 64,516 | 49 (2.14) | 18 (0.80) | 37 (1.66) | 16 (0.74) | 13 (0.60) | |

128 | 260,100 | 49 (8.93) | 18 (3.23) | 36 (6.60) | 16 (3.11) | 13 (2.53) |

It should be noted that since the results in Tables 6 and 7 are using left preconditioning, the relative residuals are actually preconditioned relative residuals. So the norms of the various preconditioners could be having an effect on the iteration counts. Since all of our numerical tests were run using the method of manufactured solutions, in Tables 8 and 9 we examine the relative errors for each method.

In Table 8 we report the relative errors for various numbers of stages for Radau IIA and Lobatto IIIC with the five preconditioners , , , , and all used as left preconditioners. The problem size is fixed at , which is the largest of the problems we considered and , where is the degree of the Lagrange polynomial in space and is the order of the corresponding IRK method. The stopping criteria is the same as in all previous tables. With left preconditioning, problems with Radau IIA and Lobatto IIIC methods and all numbers of stages yielded very similar relative errors.

IRK method | |||||
---|---|---|---|---|---|

R IIA 2 | |||||

R IIA 3 | |||||

R IIA 4 | |||||

R IIA 5 | |||||

R IIA 6 | |||||

L IIIC 2 | |||||

L IIIC 3 | |||||

L IIIC 4 |

In Table 9 we report the relative errors for various numbers of stages for Radau IIA and Lobatto IIIC with the four preconditioners , , , and this time all used as right preconditioners. We do not include the optimized GSL-based preconditioner since it is optimized for left preconditioning. The problem size is again fixed at , and as before , where is the degree of the Lagrange polynomial in space and is the order of the corresponding IRK method. The stopping criteria is the same as in all previous tables. With right reconditioning, our preconditioner achieves almost an order of magnitude better relative error on almost all of the runs.

IRK method | ||||
---|---|---|---|---|

R IIA 2 | ||||

R IIA 3 | ||||

R IIA 4 | ||||

R IIA 5 | ||||

R IIA 6 | ||||

R IIA 7 | ||||

L IIIC 2 | ||||

L IIIC 3 | ||||

L IIIC 4 | ||||

L IIIC 5 |

## 6 Conclusion

In this paper, we have introduced a new preconditioner for the large, structured systems appearing in implicit Runge–Kutta time integration of parabolic PDE. Our preconditioner is based on a block factorization, and for scalability, we have used a single AMG V-cycle on all subsolves. We have compared our preconditioner in condition number and eigenvalue distribution to other preconditioners, and have demonstrated its effectiveness on the heat equation as a simple test problem. We have found that it is scalable (-independent) and yields better timing results than other preconditioners currently in the literature: block Jacobi, block Gauss–Seidel, and the optimized block Gauss–Seidel method of [staff2006preconditioning]. We ran experiments with implicit Runge–Kutta stages up to , and have found that the new preconditioner outperforms the others, with the improvement becoming more pronounced as spatial discretization is refined and as temporal order is increased.

Comments

There are no comments yet.