Stage-parallel fully implicit Runge-Kutta solvers for discontinuous Galerkin fluid simulations Will Pazner1 and Per-Olof Persson2 7 1 1Division of Applied Mathematics, Brown University, Providence, RI, 02912 0 2Department of Mathematics, University of California, Berkeley, Berkeley, CA, 94720-3840 2 n a J Abstract 5 In this paper, we develop new techniques for solving the large, coupled linear systems that 2 arisefromfullyimplicitRunge-Kuttamethods. Thismethodmakesuseoftheiterativeprecon- ] ditionedGMRESalgorithmforsolvingthelinearsystems,whichhasseensuccessforfluidflow A problemsanddiscontinuousGalerkindiscretizations. Bytransformingtheresultinglinearsys- N temofequations,onecanobtainamethodwhichismuchlesscomputationallyexpensivethan theuntransformedformulation,andwhichcomparescompetitivelywithothertime-integration . h schemes,suchasdiagonallyimplicitRunge-Kutta(DIRK)methods. Wedevelopandtestsev- t a eral ILU-based preconditioners effective for these large systems. We additionally employ a m parallel-in-timestrategytocomputetheRunge-Kuttastagessimultaneously. Numericalexper- imentsareperformedontheNavier-StokesequationsusingEulervortexand2Dand3DNACA [ airfoil test cases in serial and in parallel settings. The fully implicit Radau IIA Runge-Kutta 1 methods compare favorably with equal-order DIRK methods in terms of accuracy, number of v GMRES iterations, number of matrix-vector multiplications, and wall-clock time, for a wide 1 range of time steps. 8 1 Keywords: implicit Runge-Kutta; discontinuous Galerkin; preconditioned GMRES; parallel-in-time 7 0 . 1 Introduction 1 0 7 ThediscontinuousGalerkinmethod,introducedin1973byReedandHillfortheneutrontransport 1 equation [29], has seen in recent years increased interest for fluid dynamics applications [24]. The : v discontinuous Galerkin method is a high-order finite element method suitable for use on unstruc- i tured meshes with polynomials of arbitrarily high degree. For many fluid flow problems, explicit X time integration methods have the downside of restrictive time step conditions, partly because of r a the use of high-degree polynomials but more fundamentally because of the need to employ highly graded and/or anisotropic elements for many realistic flow problems [27]. Therefore, for many applications it is desirable to use an implicit time integration scheme. Implicit time integration methods for DG have been much studied. Multi-step backward differ- entiation formulas (BDF) and single-step diagonally implicit Runge-Kutta (DIRK) methods have been applied to discontinuous Galerkin discretizations for fluid flow problems [26, 27]. Nigro et. al have seen success applying multi-stage, multi-step modified extended BDF (MEBDF) and two im- plicit advanced step-point (TIAS) schemes to the compressible Euler and Navier-Stokes equations 1 [21, 22]. Additionally, in [3], Bassi et al. have used linearly implicit Rosenbrock-type to integrate DGdiscretizationsforvariousfluidflowproblems. TheBDFandDIRKmethodshavesomelimita- tions: BDFschemescanbeA-stableonlyuptosecond-order(thefamoussecondDahlquistbarrier) [10], a severe limitation when used in conjunction with a high-order spatial discretization. On the other hand, there exist high-order A-stable (and even L-stable) DIRK schemes, but these methods have a low stage-order, often resulting in order reduction when applied to stiff problems [13]. The Radau IIA methods, one class of the so-called fully implicit Runge-Kutta (IRK) methods, are high-order, L-stable, and have relatively high stage order. Consequently, these methods suffer less from order reduction than the corresponding DIRK methods when applied to stiff problems. Furthermore, these methods require only a small number of stages s, with the order of accuracy given by p = 2s−1. These methods have the drawback that each step involves the solution of large, coupled linear systems of equations. The difficulty in efficiently implementing such methods has caused them to remain not widely used or studied for practical applications [8, 9]. There has been previous work on improving the efficiency of solving these large, coupled systems. In [16], Jay and Braconnier develop a parallelizable preconditioner for IRK methods by means of Hairer and Wanner’s W-transformation. In [11], De Swart et al. have developed a parallel software package for the four-stage Radau IIA method, and Burrage et at. have developed a matrix-free, parallel implementation of the fifth-order Radau IIA method in [6]. Inthis paper, wedevelop anew strategyfor efficientlysolving theresulting largelinearsystems by means of the iterative preconditioned GMRES algorithm. A simple transformation of the lin- ear system results in a significant reduction of the cost per GMRES iteration. Furthermore, the block ILU(0) preconditioner, used successfully with implicit time-integrators for the discontinuous Galerkin method in [28], proves to be effective also for these large systems. A shifted, uncoupled, block ILU(0) factorization is also found to be an effective preconditioner, with the advantage of allowing parallelism in time by computing the stage solutions simultaneously. The structure of this paper is as follows. In Section 2, we describe the governing equations and DG spatial discretization. In Section 3, we discuss the time integration schemes used in this paper. Then, in Section 4, we introduce the transformation used to reduce the solution cost, and discuss the preconditioners used for the GMRES method. Finally, in Section 5, we perform numerical experiments on a variety of test cases, in two and three spatial dimensions. 2 Equations and spatial discretization The equations considered are the time-dependent, compressible Navier-Stokes equations, ∂ρ ∂ + (ρu )=0 (1) j ∂t ∂x j ∂ ∂ ∂p ∂τ ij (ρu )+ (ρu u )+ = for i=1,2,3, (2) i i j ∂t ∂x ∂x ∂x j i j ∂ ∂ ∂q ∂ (ρE)+ (u (ρE+p))=− j + (u τ ), (3) j i ij ∂t ∂x ∂x ∂x j j j where ρ is the density, u is the ith component of the velocity, and E is the total energy. The i viscous stress tensor and heat flux are given by ∂u ∂u 2∂u µ ∂ p 1 τ =µ i + j − kδ and q =− E+ − u u , (4) ij ij j k k ∂x ∂x 3∂x Pr∂x ρ 2 (cid:18) j i k (cid:19) j (cid:18) (cid:19) 2 where µ is the viscosity coefficient, and Pr = 0.72 is the Prandtl number. For an ideal gas, the pressure p is given by the equation of state 1 p=(γ−1)ρ E− u u , (5) k k 2 (cid:18) (cid:19) where γ is the adiabatic gas constant. For the viscous problems, we introduce an isentropic as- sumption of the form p = Kργ, for a given constant K, as described in [17]. This additional simplification can be thought of as an artificial compressibility model for the incompressible flows simulatedinSections5.2and5.3. Weprefertosolvethecompressibleequationsbecausetheyresult in a system of ODEs rather than differential-algebraic equations, and thus do not need specialized projection-type solvers. This model decouples equation (3) from equations (1) and (2), and there- fore results in one fewer component to solve for. We remark that this simplification does not result in a significant difference in the relative performance of the time integrators studied in this paper, asshowninSection(5.1),wherethefullcompressibleEulerequationsaresolved,withnoisentropic assumption. We rewrite equations (1), (2), and (3) in the form ∂u +∇·F (u)−∇·F (u,∇u)=0, (6) i v ∂t where u is a vector of the conserved variables, and F ,F are the inviscid and viscous flux func- i v tions, respectively. The spatial domain Ω is discretized into a triangulation, and the solution u is approximated by piecewise polynomials of a given degree. Equation (6) is discretized by means of the discontinuous Galerkin method, where the viscous terms are treated using the compact DG (CDG) scheme [23]. Using a nodal basis function for each component, we write the global solution vector as u, and obtain a semi-discrete system of ordinary differential equations of the form ∂u M =f(u), (7) ∂t where M is the mass matrix, and f is a nonlinear function of the n unknowns u. The standard method of lines approach allows for the solution of this system of ordinary differential equations by means of a range of numerical time integrators, such as the implicit Runge-Kutta methods that are the focus of this paper. 2.1 Block structure of the Jacobian Weconsiderthevectorofunknownsutobeorderedsuchthatthemdegreesoffreedomassociated with one element of the triangulation appear consecutively. We suppose that there are a total of T elements, such that there are a total of n=Tm degrees of freedom. Then, the Jacobian matrix J can be seen as a T ×T block matrix, with blocks of size m×m. The ith row consists of blocks on the diagonal, and in columns j, where elements i and j share a common edge, such that the total number of off-diagonal blocks in the ith row is equal to the number of neighbors of element i. We note that the off-diagonal blocks of size m×m are themselves sparse, but for the sake of simplicity we will consider them as dense matrices. The mass matrix M is a T ×T block diagonal matrix, with blocks of size m×m, and therefore matrices of the form αM −βJ have the same sparsity pattern as the Jacobian. 3 3 Time integration In this paper, we will focus on the one-step, multi-stage Runge-Kutta methods. Given initial conditionsu =u(t ),agenerals-stage,pth-orderRunge-Kuttamethodforadvancingthesolution 0 0 to u =u(t +∆t)+O(∆tp+1) can be written as 1 0 s Mk =f t +∆tc ,u +∆t a k , (8) i 0 i 0 ij j   j=1 (cid:88)  s  u =u +∆t b k , (9) 1 0 i i i=1 (cid:88) where the coefficients a ,b , and c can be expressed compactly in the form of the Butcher tableau, ij i i c a ··· a 1 11 1s ... ... ... ... = c A . cs as1 ··· ass bT b ··· b 1 s If the matrix of coefficients A is strictly lower-triangular, then each stage k only depends on i the preceding stages, and the method is called an explicit Runge-Kutta method. In this case, each stage may be computed by simply evaluating the function f. If A is not strictly lower-triangular, the method is called an implicit Runge-Kutta method (IRK). A particular class of implicit Runge- Kutta methods is those for which the matrix A is lower-triangular. Such methods are called diagonally-implicit Runge-Kutta (DIRK) methods [1]. Implicit Runge-Kutta methods enjoy high accuracy and very favorable stability properties, but computing the stages requires the solution of (in general nonlinear) systems of equations. In the case of DIRK methods, since A is lower- triangular, each stage k depends only on those stages k , j ≤ i, requiring the sequential solution i j of s systems, each of size n. In contrast, general IRK methods couple all of the stages, resulting in one nonlinear system of equations of size s×n. For the solution of stiff systems of equations, we are interested in those methods that are L-stable, meaning that their stability region includes the entire left half-place (A-stability), to- gether with the additional criterion that the stability function R(z) satisfies lim R(z) = 0. In z →∞ the present study, we compare the efficiency and effectiveness of several L-stable IRK and DIRK schemes. ThemethodsconsideredarelistedinTable1. TheIRKschemesconsideredaretheRadau IIA schemes, which are L-stable, s-stage schemes of order 2s−1 based on the Radau right quadra- ture. Theconstructionoftheseschemescanbefoundin[15]. Thetwo-stageandthree-stageRadau IIA methods are listed as RADAU23 and RADAU35, respectively. The DIRK schemes considered are the three-stage, third-order L-stable scheme denoted DIRK33, which is derived in detail in [1], and the six-stage, fifth-order scheme constructed in [5], and denoted ESDIRK65. The latter scheme is an explicit singly diagonally implicit Runge-Kutta (ESDIRK) method, meaning that the first diagonal entry of the Butcher matrix is zero, and the remaining diagonally entries are nonzero and equal. In addition to the third- and fifth-order methods, we also consider the seventh- and ninth-order Radau IIA methods, although we do not compare these methods to equal-order DIRK methods. The Butcher tableaux for the methods considered are given in Appendix A. Also of interest is the phenomenon of order reduction, whereby, when applied to stiff problems, the overall order of accuracy is reduced from p to the stage order of the method (denoted q) [13]. 4 Table 1: L-stable implicit Runge-Kutta schemes considered Scheme Order Total stages Implicit stages Stage order Leading error coefficient RADAU23 3 2 2 2 1.39×10 2 − DIRK33 3 3 3 1 2.59×10 2 − RADAU35 5 3 3 3 1.39×10 4 − ESDIRK65 5 6 5 2 5.30×10 4 − RADAU47 7 4 4 4 7.09×10 7 − RADAU59 9 5 5 5 2.19×10 9 − The stage order q is defined as q =min{p,q }, for i=1,...,s, where q is defined by i i s u(t +∆tc )=u +∆t a k +O(∆tqi+1). (10) 0 i 0 ij j j=1 (cid:88) It can be shown that the maximum stage order for any DIRK method is 2 [14], whereas the stage order for the Radau IIA methods is given by the number of stages, q = s [19]. The DIRK33 has stage order of q =1. An advantage of the ESDIRK methods such as ESDIRK65 is that they have stage order of q =2 [18]. TheRadauIIAmethodsareveryattractivebecauseoftheirhighorderofaccuracy,smallnumber of stages, high stage order, and L-stability, but solving the coupled system of s×n equations is computationallyexpensive. Supposingthatwesolvethenonlinearsystemofequationsforthestages k bymeansofNewton’smethod, thenateachiterationwemustsolvealinearsystemofequations i byinvertingtheJacobianmatrixoftheright-handside,f(t,u). AssumingadenseJacobianmatrix, and solution via Gaussian elimination (or LU factorization), then the cost of performing a linear solvescalesasthecubeofthenumberofunknowns. Therefore,thecostperlinearsolveforageneral IRK method is O(s3n3), whereas the cost per solve for a DIRK method scales like O(sn3). In [7], Butcherdescribeshowtotransformtheresultingsetoflinearequationstoreducethecomputational work for solving the IRK systems to O(2sn3). Despite this reduction in computational complexity, the cost of solving the large systems of equations has proven in practice to be prohibitive [8]. On the other hand, DIRK methods have proven be popular and effective for solving computational fluid dynamics problems [4], at the cost of lower stage order and an increased number of stages. 4 Efficient solution of implicit Runge-Kutta systems In this section we describe a method for efficiently solving the systems arising from general IRK methods when applied to discontinuous Galerkin discretizations. The Jacobian matrices of the function f are sparse, block-structured matrices, which lend themselves to solution via iterative Krylov subspace methods. In particular, we consider the solution of these systems by means of the GMRES method with a zero fill-in block ILU(0) preconditioner, as in [28]. In this case, each iteration of the GMRES method requires one matrix-vector multiplication, and one application of the ILU(0) preconditioner. In order to efficiently solve the linear systems resulting from IRK methods, we will rewrite the system of equations in such a way so as to reduce the cost of a matrix-vector multiplication from s2n2 to order sn2. 5 4.1 Transformation of the system of equations Recalling equation that the stages k are given by the equation i s Mk =f t +∆tc ,u +∆t a k , (11) i 0 i 0 ij j   j=1 (cid:88)   we define K to be the concatenation of the vectors k , U to be the concatenation of s copies of i 0 u ,andF thefunctionf appliedcomponent-wiseonthesevectors. Then,werewriteequation(11) 0 in vector form as (I ⊗M)K =F (t +∆tc,U +∆t(A⊗I )K), (12) s 0 0 n where ⊗ is the Kronecker product, and I is the n×n identity matrix. n This nonlinear system can be solved by means of Newton’s method, which will require solving at each step a linear system of the form M 0 a J ··· a J k r 11 1 1s 1 1 1  ... −∆t ... ... ...  ... = ...  (13) 0 M a J ··· a J , k r    s1 s ss s  s   s         for the residual vectors (r ,...,r )T, where the matrices on the left-hand side are s × s block 1 s matrices blocks, with each block of size n×n. We use the following notation for the Jacobian matrix of f, s J =J t +∆tc ,u +∆t a k . i f 0 i 0 ij j   j=1 (cid:88)   We can rewrite equation (13) in the following form, a J ··· a J 11 1 1s 1 Is⊗M −∆t ... ... ... K =R, (14) a J ··· a J   s1 s ss s     The sparsity pattern of the matrix in (14) is simply that of the Jacobian matrix J repeated s×s times, and we can conclude that the cost of computing a matrix-vector product with this matrix is s2 times that of computing the matrix-vector product of one Jacobian matrix. In order to reduce the cost of the matrix-vector multiplication, we perform a simple transfor- mation in to rewrite (14) in a slightly modified form. We begin by defining s w = a k , i ij j j=1 (cid:88) and similarly, we let W denote the vectors w stacked, such that i W =(A⊗I )K. n Then, we rewrite the nonlinear system of equations (11) in terms of the variables w to obtain i Mk =f(t +∆tc ,u +∆tw ), i 0 i 0 i 6 or, equivalently, in the case that the Butcher matrix A is invertible, (A 1⊗M)W =F(t +∆tc,U +∆tW). (15) − 0 0 In the transformed variables, the new solution u can be written as 1 u =u +∆t(bTA 1⊗I )W, 1 0 − n In the case of the Radau IIA methods, bTA 1 =(0,...,0,1), and so this further simplifies to − u =u +∆tw . 1 0 s Solving equation (15) with Newton’s method gives rise to the linear system of equations J 0 ··· 0 1 0 J ··· 0 2 A−1⊗M −∆t ... ... ... ... W =R. (16)      0 0 ··· Js        Theadvantageofthisformulationisthattheresultingsystemenjoysgreatersparsity. Theresulting matrixisas×sblockmatrix,withmultiplesofthemassmatrixineveryoff-diagonalblock,andwith matrices of the form (A 1) M −∆tJ along the diagonal. Computing a matrix-vector product − ii i with this s×s block matrix requires performing s matrix-vector multiplications with the mass matrix, and s matrix-vector multiplications with a Jacobian J . Therefore the cost of computing i such products scales as s times the cost of computing one matrix-vector product with the Jacobian matrix. We additionally remark that the fully-implicit IRK methods requiring storing each of the s Jacobian matrices J , resulting in memory usage that is s-times that of the DIRK methods. A i further advantage of the transformed system of equations is that the memory requirements for the ILU-based preconditioners are reduced, as discussed in the following sections. 4.2 Preconditioning The use of an appropriate preconditioner is essential in accelerating the convergence of a Krylov subspace method such as GMRES. We briefly describe the block ILU(0) preconditioner from [28]. 4.2.1 Block ILU(0) preconditioner The block ILU(0) (or zero fill-in) preconditioner is a method for obtaining block lower- and upper- triangular matrices L˜ and U˜ given a block sparse matrix B. These matrices are obtaining by computing the standard block LU factorization, but discarding any blocks which do not appear in the sparsity pattern of B. We denote the block in the ith row and jth column as B . The ILU(0) ij algorithm can be written as shown in Algorithm 1. If we impose the condition on the triangulation of our domain that, if elements j and k both neighborelementi,thenelementsjandkarenotneighborsofeachother,thentheILU(0)algorithm hastheparticularlysimpleform,showinAlgorithm2. Inpractice,mostwell-shapedmeshessatisfy this condition and henceforth we will use this simpler algorithm. 7 Algorithm 1 Block ILU(0) algorithm 1: for i=1 to T do 2: for neighbors j of i with j >i do 3: Bji ←BjiBi−i1 4: Bjj ←Bjj −BjiBij 5: for neighbors k of j and i with k >j do 6: Bjk ←Bjk−BjiBik 7: L˜ ←I+strict block lower triangle of B 8: U˜ ←block upper triangle of B Algorithm 2 Simplified block ILU(0) algorithm 1: for i=1 to T do 2: for neighbors j of i with j >i do 3: Bji ←BjiBi−i1 4: Bjj ←Bjj −BjiBij 5: L˜ ←I+strict block lower triangle of B 6: U˜ ←block upper triangle of B 4.3 Preconditioning the large system In the case of the general IRK methods, we must solve systems of the form J 0 ··· 0 1 0 J ··· 0 2 BW =R, B =A−1⊗M −∆t ... ... ... ... . (17)      0 0 ··· Js        We remark that the matrix B can now be considered as a s×s block matrix, with blocks of size n×n. Each n×n block is itself a T ×T block matrix, with subblocks of size m×m. We introduce the notation B to denote the (i,j) subblock of the (k,(cid:96)) block of B. That is to say, B is k(cid:96),ij k(cid:96),ij the (i,j) block of the matrix (A)−k(cid:96)1M −δk(cid:96)∆tJk, where δk(cid:96) is the Kronecker delta. 4.3.1 Stage-coupled block ILU(0) preconditioner We consider two preconditioners for this large sn×sn system. The first is the standard (stage- coupled)blockILU(0)preconditioner,whichcanbecomputedusingAlgorithm3. Wenotethatthis preconditioner couples all s stages of the method. This preconditioner requires storing s Jacobian- sized diagonal blocks, and s2−s off-diagonal blocks, each the same size as the mass matrix. 8 Algorithm 3 ILU(0) algorithm for IRK systems of the form (17) 1: for k =1 to s do 2: for i=1 to T do 3: for neighbors j of i with j >i do 4: Bkk,ji ←Bkk,jiBk−k1,ii 5: Bkk,jj ←Bkk,jj −Bkk,jiBkk,ij 6: for (cid:96)=k+1 to s do 7: B(cid:96)k,ii ←B(cid:96)k,iiBk−k1,ii 8: B(cid:96)(cid:96),ii ←B(cid:96)(cid:96),ii−B(cid:96)k,iiBkk,ii 9: L˜ ←I+strict block lower triangle of B 10: U˜ ←block upper triangle of B 4.3.2 Stage-uncoupled, shifted ILU(0) preconditioner In order to avoid the above coupling of the stages, we can compute a simplified preconditioner in the form of the following block matrix, L˜ U˜ 0 ··· 0 1 1 0 L˜ U˜ ··· 0 2 2  ... ... ... ... , (18)    0 0 ··· L˜ U˜   s s    whereL˜iU˜i istheblockILU(0)factorizationofamatrixoftheform A−ii1+αi M−∆tJi. Welet α denote a shift, so that the standard unshifted factorization corresponds to α =0. The so-called i i (cid:0) (cid:1) shifted ILU factorization, described by Manteuffel in [20], can result in eigenvalues clustered away from the origin, and hence faster convergence in GMRES. Indeed, our experience shows that the unshifted preconditioner underperforms certain other choices of shift. The preconditioner has several advantages over the stage-coupled block ILU(0) preconditioner. Thefirstisthatitiseasilyconstructedusinganalready-implementedblockILU(0)factorizationof the Jacobian matrix. The second is that none of the stages are coupled, allowing for both efficient computation and application. In particular, this has implications for the parallelization of the preconditioner, which we discuss in Section 4.5. Finally, as this preconditioner does not include any off diagonal blocks, the memory requirements are exactly s times that of the standard block ILU(0) used for the DIRK methods. As mentioned, the unshifted preconditioner, with α =0 for all i, such that L˜ U˜ is the ILU(0) i i i factorizationoftheithdiagonalblockofthematrixB,isanaturalchoice. Thischoiceofcoefficients ignoresalltheoff-diagonalmassmatrices. Bymakingcertainjudiciouschoicesofthecoefficientsα , i we can attempt to compensate for the off-diagonal blocks by adding multiples of the mass matrix back to the diagonal entries. In particular, our numerical experiments have shown that setting αi = j=i A−ji1 results in a more efficient preconditioner, requiring fewer GMRES iterations in order to c(cid:54) onverge to a given desired tolerance. (cid:80) (cid:12) (cid:12) (cid:12) (cid:12) 9 4.4 Computational cost and memory requirements InordertocomparethecomputationalcostofthetransformedIRKimplementationdescribedabove withboththatoftheuntransformedformulation,andwiththeusualDIRKmethods,wesummarize the computational cost associated with solving the resulting linear systems. We note that the IRK methods require the solution of one large, coupled system, whereas the DIRK methods require the solution of s smaller systems. In Table 2 we record the leading terms of the computational cost of operationsrequiredtobeperformedeveryiteration. WerecallthatsisthenumberofRunge-Kutta stages, m is the number of degrees of freedom per mesh element, T is the total number of elements in the mesh, and r is the number of neighbors per element. Here we assume that the m×m blocks oftheJacobianmatrixaredense,andhencerequire2m2 floatingpointoperationspermatrix-vector multiply. ComputingthepreconditionerrequirestheLUfactorizationofthediagonalblocks,which incursacostof 2m3floatingpointoperationsperblock. EachiterationinNewton’smethodrequires 3 re-evaluation of the Jacobian matrix, and therefore also the re-computation of the preconditioner. Hence, thepreconditionermustbecomputedonceperlinearsolve. Thecostsassociatedwiththese calculations are listed in Table 3. Table 2: Per-iteration computational costs for solving implicit Runge-Kutta systems Operation Cost (leading term) Untransformed IRK matrix-vector product s2m2(r+1)T Transformed IRK matrix-vector product sm2(r+s)T DIRK matrix-vector product m2(r+1)T Coupled preconditioner application (IRK) sm2(r+s)T Uncoupled preconditioner application (IRK) sm2(r+1)T Preconditioner application (DIRK) m2(r+1)T Table 3: Computational cost of computing the preconditioner (once per solve) systems Operation Cost (leading term) Computing coupled block ILU(0) preconditioner (IRK) s(m3+(r+s)m2)T Computing uncoupled ILU(0) preconditioner (IRK) s(m3+rm2)T Computing block ILU(0) preconditioner (DIRK) (m3+rm2)T We remark that each GMRES iteration using the formulation described in Section 4 requires a factor of s fewer floating-point operations per iteration than the na¨ıve IRK implementation. The stage-uncoupled IRK preconditioner is less expensive to both compute and apply than the stage- coupledpreconditioner. Wealsonotethatforequalorderofaccuracy,theRadauIIAIRKmethods require fewer implicit stages than do the DIRK methods. Each such implicit stage incurs the cost of assembling the Jacobian matrix. This cost is problem-dependent, but is in general non-trivial, and in the model problems considered in this paper, it scales like O(m3T). Finally, we present the memory requirements for the IRK and DIRK methods, and the stage- coupledanduncoupledpreconditionersinTable4. WenotethatforthetransformedIRKmethods, 10