A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE

10/22/2020 ∙ by Md Masud Rana, et al. ∙ Texas Tech University 0

A new preconditioner based on a block LDU factorization with algebraic multigrid subsolves for scalability is introduced for the large, structured systems appearing in implicit Runge-Kutta time integration of parabolic partial differential equations. This preconditioner is compared in condition number and eigenvalue distribution, and in numerical experiments with others in the literature: block Jacobi, block Gauss-Seidel, and the optimized block Gauss-Seidel method of 10.4173/mic.2006.2.3. Experiments are run with implicit Runge-Kutta stages up to s=7, and it is 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.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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
Table 1: Common IRK Methods

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 form

In 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
Table 2: Condition numbers of right preconditioned matrices with preconditioners , , , and applied to a 2D problem with -stage Radau IIA methods. Here, and , where is the degree of the Lagrange polynomial basis functions in space. Preconditioners are constructed exactly.

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.

Figure 1: Eigenvalues of the matrices , , , , and for 2D problem with Radau IIA (a), (b), (c), and (d). The -axis is the real part and the -axis is the imaginary part of the eigenvalue. Preconditioners are constructed exactly.

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
Table 3: Condition numbers of the left-preconditioned system with preconditioners ( with optimal coefficients) and for various IRK methods applied to and problems. Here, and , where is the degree of the Lagrange polynomial in space and is the order of the corresponding IRK method. Preconditioners are constructed exactly.

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 .

Figure 2: Eigenvalues of the matrix , , and for 2D problem with Radau IIA (a), (b), (c), and (d). The -axis is the real part and the -axis is the imaginary part of the eigenvalue. Preconditioners are constructed exactly.
Figure 3: Eigenvalues of the matrix , , and for 2D problem with Lobatto IIIC stages (a) and (b). The -axis is the real part and the -axis is the imaginary part of the eigenvalue. Preconditioners are constructed exactly.

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)
Table 4: Iteration counts and elapsed time (times in seconds are shown in parentheses) for right-preconditioned GMRES to converge with relative tolerance for a problem with -stage Radau IIA methods with preconditioners , , , and . Here we choose , where is the degree of the Lagrange polynomial in space. Preconditioners are approximated using one AMG V-cycle for each subsolve.

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)
Table 5: Iteration counts and elapsed time (times in seconds are shown in parentheses) for right-preconditioned GMRES to converge with relative tolerance for a problem with -stage Lobatto IIIC methods with preconditioners , , , and . Here we choose , where is the degree of the Lagrange polynomial in space. Preconditioners are approximated using one AMG V-cycle for each subsolve.

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)
Table 6: Iteration counts and elapsed time (times in seconds are shown in parentheses) for left-preconditioned GMRES to converge with relative tolerance for a problem with -stage Radau IIA methods with preconditioners , , , , and . Here we choose , where is the degree of the Lagrange polynomial in space. Preconditioners are approximated using one AMG V-cycle for each subsolve.

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)
Table 7: Iteration counts and elapsed time (times in seconds are shown in parentheses) for left-preconditioned GMRES to converge with relative tolerance for a problem with -stage Lobatto IIIC methods with preconditioners , , , , and . Here we choose , where is the degree of the Lagrange polynomial in space. Preconditioners are approximated using one AMG V-cycle for each subsolve.

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
Table 8: Relative error for left-preconditioned GMRES converging to a relative tolerance of for a problem with various IRK methods with preconditioners , , , , and . Here, and , where is the degree of the Lagrange polynomial in space and is the order of the corresponding IRK method. Preconditioners are approximated using one AMG V-cycle for each subsolve.

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
Table 9: Relative error for right-preconditioned GMRES converging to a relative tolerance of for a problem with various IRK methods with preconditioners , , , and . Here, and , where is the degree of the Lagrange polynomial in space and is the order of the corresponding IRK method. Preconditioners are approximated using one AMG V-cycle for each subsolve.

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.

References