On the convergence of iterative solvers for polygonal discontinuous Galerkin discretizations
OON THE CONVERGENCE OF ITERATIVE SOLVERS FOR POLYGONALDISCONTINUOUS GALERKIN DISCRETIZATIONS
WILL PAZNER, PER-OLOF PERSSON
Abstract.
We study the convergence of iterative linear solvers for discontinuous Galerkin discretizations ofsystems of hyperbolic conservation laws with polygonal mesh elements compared with that of traditional triangularelements. We solve the semi-discrete system of equations by means of an implicit time discretization method, usingiterative solvers such as the block Jacobi method and GMRES. We perform a von Neumann analysis to analyticallystudy the convergence of the block Jacobi method for the two-dimensional advection equation on four classes of regularmeshes: hexagonal, square, equilateral-triangular, and right-triangular. We find that hexagonal and square meshesgive rise to smaller eigenvalues, and thus result in faster convergence of Jacobi’s method. We perform numericalexperiments with variable velocity fields, irregular, unstructured meshes, and the Euler equations of gas dynamics toconfirm and extend these results. We additionally study the effect of polygonal meshes on the performance of blockILU(0) and Jacobi preconditioners for the GMRES method.
1. Introduction.
In recent years, the Discontinuous Galerkin (DG) method has become apopular choice for the discretization of a wide range of partial differential equations [27, 6, 15].This is partly because of its many attractive properties, such as the arbitrarily high degrees ofapproximation, the rigorous theoretical foundation, and the ability to use fully unstructured meshes.Also, due to its natural stabilization mechanism based on approximate Riemann solvers, it has inparticular become widely used in fluid dynamics applications where the high-order accuracy isbelieved to produce improved accuracy for many problems [32].Most work on DG methods has been based on meshes of either simplex elements (trianglesand tetrahedra), block elements (quadrilaterals and hexahedra), or combinations of these such asprism elements. This is likely because of the availability of excellent automatic unstructured meshgenerators, at least for the simplex case [22, 28, 30], and also because of the advantages with theouter-product structure of block elements. However, it is well known that since no continuity is en-forced between the elements, it is straightforward to apply the DG methods to meshes with elementsof any shapes (even non-conforming ones). For example, vertex-centered DG methods based on thepolygonal dual meshes were studied in [5, 18]. This is a major advantage over standard continuousFEM methods, which need significant developments for the extension to arbitrary polygonal andpolyhedral elements [19].In the finite volume CFD community, there has recently been considerable interest in meshesof arbitrary polygonal and polyhedral elements. In fact, the popular vertex-centered finite volumemethod applied to a tetrahedral mesh can be seen as a cell-centered method on the dual polyhedralmesh. Because of this, a number of methods have been proposed for generation of polyhedralmeshes, which in many cases have advantages over traditional simplex meshes [21, 12]. Althoughit is still unclear exactly what benefits these elements provide, they have been reported to be bothmore accurate per degree of freedom and to have better convergence properties in the numericalsolvers than for a corresponding tetrahedral mesh [23, 2]. There have also been studies showingthat vertex-centered schemes are preferred over cell-centered [10, 9], again indicating the benefitsof polyhedral elements.Inspired by the promising results for polyhedral finite volume method, and the fact that DG isa natural higher-order extension of these schemes, in this work we study some of the properties ofDG discretizations on polygonal meshes. To limit the scope, we only investigate the convergenceproperties of iterative solvers for the discrete systems, assuming an equal number of degrees offreedom per unit area for all element shapes. Future work will also investigate the accuracy of the a r X i v : . [ m a t h . NA ] J a n olutions on the different meshes. We first consider the iterative block-Jacobi method applied toa pure convection problem, which in the constant coefficient case can be solved analytically usingvon Neumann analysis. Next we apply the solver to Euler’s equations of gas dynamics for relevantmodel flow problems, to obtain numerical results for the convergence of the various element shapes.We consider regular meshes of hexagons, squares, and two different configurations of triangles, aswell as the dual of fully unstructured triangular Delaunay refinement meshes. We also performnumerical experiments with the GMRES Krylov subspace solver and a block-ILU preconditioner.Although the results are not entirely conclusive, most of the results indicate a clear benefit withthe hexagonal and quadrilateral elements over the triangular ones.The paper is organized as follows. In Section 2, we describe the spatial and the temporaldiscretizations, and introduce the iterative solvers. In Section 3 we perform the von Neumannanalysis of the constant coefficient advection problem, in 1D and for several mesh configurationsin 2D. In Section 4 we show numerical results for more general advection fields, for more generalmeshes, as well as for the Euler equations and the GMRES solver. We conclude with a summaryof our findings as well as directions for future work.
2. Numerical methods.2.1. The discontinuous Galerkin formulation.
We consider a system of m hyperbolicconservation laws given by the equation(1) (cid:40) ∂ t u + ∇ · F ( u ) = 0 , ( t, x ) ∈ [0 , T ] × Ω u (0 , x ) = u ( x ) . In order to describe the discontinuous Galerkin spatial discretization, we divide the spatial domainΩ ⊆ R into a collection of elements, to form the triangulation T h = { K i } . Often the elements K i are considered to be triangles or quadrilaterals, but in this paper we allow the elements to bearbitrary polygons in order to study the impact of different tessellations on the efficiency of thealgorithm.Let V h = (cid:110) v h ∈ L (Ω) : v h (cid:12)(cid:12) K i ∈ P p ( K i ) (cid:111) denote the space of piecewise polynomials of degree p . We let V mh denote the space of vector-valued functions of length m , with each component in V h .Note that continuity is not enforced between the elements. We derive the discontinuous Galerkinmethod by replacing u in equation (1) by an approximate solution u h ∈ V mh , and then multiplyingequation by a test function v h ∈ V mh . We then integrate by parts over each element. Becausethe approximate solution u h is potentially discontinuous at the boundary of an element, the fluxfunction F is approximated by a numerical flux function (cid:98) F , which takes as arguments u + , u − ,and n , denoting the solution on the exterior and interior of the element, and the outward-pointingnormal vector, respectively. Then, the discontinuous Galerkin method reads:Find u h ∈ V mh such that, for all v h ∈ V mh , (cid:90) K i ∂ t u h · v h dx − (cid:90) K i F ( u h ) : ∇ v h dx + (cid:73) ∂K i (cid:98) F ( u + , u − , n ) · v h ds = 0 . (2) As a first example, we consider the two-dimensional scalar advec-tion equation(3) u t + ∇ · ( β u ) = 0 , or a given (constant) velocity vector β = ( α, β ). We solve this equation in the domain [0 , π ] × [0 , π ], with periodic boundary conditions. The exact solution to this equation is given by(4) u ( t, x, y ) = u ( x − αt, y − βt ) , where u is the given initial state.In order to define the discontinuous Galerkin method for equation (3), we define the upwindflux by(5) (cid:98) F ( u + , u − , n ) = (cid:40) u − if β · n ≥ u + if β · n < u h as a vector U consisting of the coefficients ofthe expansion of u h in terms of an orthogonal Legendre polynomial modal basis of the functionspace V mh . Discretizing equation (3) results in a linear system of equations, which we can write as(6) M ( ∂ t U ) + L U = 0 , where the mass matrix M corresponds to the first term on the left-hand side of (2), and L consistsof the second two terms on the left-hand side. The mass matrix is block-diagonal, and the matrix L is a block matrix, with blocks along the diagonal, and off-diagonal blocks corresponding to theboundary terms from the neighboring elements. We consider the solution of (6) by meansof implicit time integration schemes, the simplest of which is the standard backward Euler scheme,(7) ( M + k L ) U n +1 = M U n . Furthermore, each stage of a higher-order scheme, such as a diagonally-implicit Runge-Kutta(DIRK) scheme [1], can be written as a similar equation. The block sparse system can be solvedefficiently by means of an iterative linear solver. In this paper, we consider two solvers: the simpleblock Jacobi method, and the preconditioned GMRES method.
A popular and simple iterative solver is the block Jacobimethod, defined as follows. Each iteration of the method for solving the linear system A x = b isgiven by(8) x ( n +1) = D − b + R J x ( n ) , where D is the block-diagonal part of A , and R J = I − D − A . This simple method has the advan-tage that it is possible to analyze the convergence properties of the method simply by examiningthe eigenvalues of the matrix R J . An upper bound of 1 for the absolute value of the eigenvalues ofthe matrix R J is a necessary and sufficient condition in order for Jacobi’s method to converge (forany choice of initial vector x (0) ). The spectral radius of R J determines the speed of convergence. Another popular and oftentimes more efficient[3] method for solving large, sparse linear systems is the GMRES (generalized minimal residual)method [25]. As with most Krylov subspace methods, the choice of preconditioner has great impacton the efficiency of the solver [26]. A simple and popular choice of preconditioner is the block Jacobipreconditioner. Each application of this preconditioner is performed by multiplying by the inverse f the block-diagonal part of the matrix. Another, often more effective choice of preconditioner isthe block ILU(0) preconditioner [8]. This preconditioner produces an approximate block-wise LUfactorization, whose sparsity pattern is enforced to be the same as that of the original matrix. Thisfactorization can be performed in-place, and requires no more storage that the original matrix.Unlike the block Jacobi method, the block ILU(0) preconditioner can be highly sensitive to theordering of the mesh elements [11, 4]. Because of this property, it is common to combine the use ofILU preconditioners with certain orderings of the mesh elements designed to increase efficiency, suchas reverse Cuthill-McKee [7], minimum degree [20], nested dissection [13], or minimum discardedfill [25].In this paper, we focus our study on the block Jacobi method, which is simpler and moreamenable to analysis. We then perform numerical experiments using both the block Jacobi methodand the preconditioned GMRES method using ILU(0) and block Jacobi preconditioning.
3. Jacobi Analysis.
We compare tessellations of the plane by four sets of generating patterns ,each consisting of one or more polygons. We consider tessellations consisting of squares, regularhexagons, two right triangles, and two equilateral triangles. The generating patterns considered areshown in Figure 1. Each generating pattern G j consists of one or two elements, labeled K j and (cid:102) K j .We will refer to these generating patterns as S, H, R, and E for squares, hexagons, right triangles,and equilateral triangles, respectively.We are interested in computing the spectral radius of the Jacobi matrix R J that arises from thediscontinuous Galerkin discretization on the mesh resulting from tessellating the plane by each ofthe four generating patterns. For the sake of comparison, we choose the elements from each of thegenerating patters to have the same area. Therefore, if the side length of the equilateral triangle is h E = h , then the two equal sides of the isosceles right triangle have side length h R = √ √ h E , thehexagon has side length h H = √ h E , and the square has side length h S = √ h E . Then, the globalsystem will have the same number of degrees of freedom regardless of choice of generating pattern. First, we compare the efficiency of each of the four types ofgenerating patterns when used to solve the advection equation (3) with the discontinuous Galerkinspatial discretization and implicit time integration. We compute the spectral radius of the matrix R J using the classical von Neumann analysis for each of the generating patterns, in a mannersimilar to [16].Let U denote the solution vector, and let its j th component, U j , which is itself a vector, denotethe degrees of freedom in G j , the j th generating pattern. We remark that in the case of squaresand hexagons, this corresponds exactly to the degrees of freedom in the element K j , but in thecase of the triangular generating patters, this corresponds to the degrees of freedom from both ofthe elements K j and (cid:102) K j . In order to determine the eigenvalues of R J , we consider the planar wavewith wavenumber ( n x , n y ) defined by(9) U j = e i ( n x x j + n y y j ) (cid:98) U , where ( x j , y j ) are fixed coordinates in G j . Then, we let (cid:96) index the generating patterns neighboring G j , and we let δ (cid:96) = ( δ x(cid:96) , δ y(cid:96) ) = ( x j − x (cid:96) , y j − y (cid:96) ) be the offsets satisfying G j + δ (cid:96) = G (cid:96) . We canthen write the solution in each of the neighboring generating patterns as(10) U (cid:96) = e i ( n x δ x(cid:96) + n y δ y(cid:96) ) U j . j (a) Square Cartesian grid K j (b) Regular hexagons K j (cid:102) K j (c) Isosceles right triangles K j (cid:102) K j (d) Equilateral triangles Fig. 1: Examples of generating patterns G j shown with bolded lines. Neighboring elements areshown unbolded.In this case we write the semi-discrete equations (6) in the following compact form(11) M j ( ∂ t U j ) + (cid:88) (cid:96) e i ( n x δ x(cid:96) + n y δ y(cid:96) ) L j(cid:96) U j = 0 , where the summation over (cid:96) ranges over all neighboring generating patterns, M j denotes the diag-onal block of M corresponding to the j th generating pattern, and L j(cid:96) denotes the block of L in the j th row and (cid:96) th column. We can write(12) (cid:98) L j = (cid:88) (cid:96) e i ( n x δ x(cid:96) + n y δ y(cid:96) ) L j(cid:96) to further simplify and obtain(13) M j ( ∂ t (cid:98) U ) + (cid:98) L j (cid:98) U = 0 . n order to solve equation (13) using an implicit method, we consider the backward Euler-typeequation(14) ( M j + k (cid:98) L j ) (cid:98) U n +1 = M j (cid:98) U n . The Jacobi iteration matrix R J can then be written as(15) (cid:99) R J j = I − D − ( M j + k (cid:98) L j ) , where the matrix D = M j + k L jj consists of the j th diagonal block of M + k L . The eigenvaluesof the matrix (cid:99) R J j control the speed of convergence of Jacobi’s method. In the simple cases ofpiecewise constant functions ( p = 0), or in the case of a one-dimensional domain, the eigenvaluescan be computed explicitly. In the more complicated case of p ≥ To illustrate the von Neumann analysis, we consider the one-dimensionalscalar advection equation(16) u t + u x = 0on the interval [0 , π ] with periodic boundary conditions. We divide the domain into N subintervals K j , each of length h . Let U denote the solution vector, and let U j denote the degrees of freedomfor the j th interval K j . For example, if piecewise constants are used, the method is identical to theupwind finite volume method, and each U j represents the average of the solution over the interval.If piecewise polynomials of degree p are used, each U j is a vector of length p + 1.For the purposes of illustration, we choose p = 1, and let U j = ( u j, , u j, ) represent the value ofthe solution at the left and right endpoints of the interval K j . Then, the local basis on the interval K j consists of the functions(17) φ j, ( x ) = j − x/h, φ j, ( x ) = x/h − j + 1 . We remark that the upwind flux in this case is always equal to the value of the function immediatelyto the left of the boundary point:(18) (cid:104) (cid:98) F ( u + , u − , x ) v ( x ) (cid:105) jh ( j − h = u j, v j, − u j − , v j, . The entries of the j th block of the mass matrix M are given by(19) ( M j ) i(cid:96) = (cid:90) jh ( j − h φ j,i ( x ) φ j,(cid:96) ( x ) dx. Additionally, we remark that the diagonal blocks of L consist of the volume integrals and rightboundary terms given by(20) ( L jj ) i(cid:96) = φ j,i ( jh ) φ i,(cid:96) ( jh ) − (cid:90) jh ( j − h φ (cid:48) j,i ( x ) φ j,(cid:96) ( x ) dx. We let A denote the backward Euler-type operator defined by(21) A = M + k L , nd, solving the equation A x = b by means of Jacobi iterations, we define the Jacobi matrix R J by(22) R J = I − D − A , where D is the matrix consisting of the diagonal blocks of A . The entries of the diagonal blocks M j and L jj can be computed explicitly using (17) to obtain(23) M j = h h h h , L jj =
12 12 −
12 12 , D j = h + k h + k h − k h + k . In order to perform the von Neumann analysis, we seek solutions of the form U j = e inhj (cid:98) U ,which allows us to explicitly compute the form of the matrix (cid:98) L j . Recalling the compact form from(13), we obtain(24) (cid:98) L j =
12 12 − e − ihn −
12 12 . Then, the Jacobi matrix (cid:99) R J j is given by(25) (cid:99) R J j = e − ihn k (2 h +3 k ) h +4 kh +6 k − e − ihn ( h − k ) kh +4 kh +6 k , whose eigenvalues λ and λ are given by(26) λ = 0 , λ = 2 k (3 k − h ) e − ihn h + 4 hk + 6 k . Therefore, each wavenumber n from 0 to 2 π/h corresponds to an eigenvalue of the Jacobi matrix R J , and the magnitude of these eigenvalues determine the speed of convergence of Jacobi’s method.In this case, the expression(27) λ max = 2 k | h − k | h + 4 hk + 6 k determines the speed of convergence of Jacobi’s method. This expression can easily be seen to bebounded above by 1 for all positive values of h and k , therefore indicating that Jacobi’s method isguaranteed to converge, unconditionally, regardless of spatial resolution or timestep. We now turn to the analysis of the four generating patterns shown inFigure 1. The analysis proceeds along the same lines as in the one-dimensional example fromSection 3.2. As an example, we present the case of piecewise constants, for which it is possible toexplicitly compute the eigenvalues of the Jacobi matrix R J . In this case the discontinuous Galerkinformulation simplifies to the upwind finite volume method(28) (cid:90) K j ∂ t u h dx + (cid:73) ∂K j (cid:98) F ( u + , u − , n ) ds = 0 . or the sake of concreteness, we assume without loss of generality that the velocity vector β = ( α, β )satisfies α, β ≥
0. In order to explicitly write the upwind flux on the meshes consisting of hexagonsand equilateral triangles, we further assume that √ α − β ≥
0, and on the mesh consisting of righttriangles we assume that α − β ≥
0. In the case of the square and hexagonal meshes, there is onlyone degree of freedom per generating pattern, and we will write u j to represent the average value ofthe solution over the generating pattern G j . We then consider the planar wave with wavenumber( n x , n y ) given by u j = e i ( n x x j + n y y j ) (cid:98) u . In the case of the square mesh with side length h S = √ h E ,the method can be written as(29) h S ( ∂ t (cid:98) u ) = − h S (cid:0) α (1 − e − in x h S ) + β (1 − e − in y h S ) (cid:1) (cid:98) u. In this case, the mass matrix M is a diagonal matrix with h S along the diagonal, and the diagonalentries of the matrix L are given by h S ( α + β ). Therefore, the eigenvalues of the Jacobi matrix R SJ = I − D − ( M + k L ) are given by(30) λ ( R SJ ) = 1 − h S + h S k ( α + β ) (cid:0) h S + h S k (cid:0) α (1 − e − in x h S ) + β (1 − e − in y h S ) (cid:1)(cid:1) = k (cid:0) αe − in x h S + βe − in y h S (cid:1) h S + k ( α + β ) . In the case of the hexagonal mesh with side length h H = √ h E , the method is(31) 3 √ h H ( ∂ t (cid:98) u ) = − h H (cid:32) (cid:16) √ α + β (cid:17) + (cid:16) − √ α + β (cid:17) e ih H (cid:16) − n x + √ n y (cid:17) + (cid:16) − √ α − β (cid:17) e ih H (cid:16) − n x − √ n y (cid:17) − βe − ih H √ n y (cid:33)(cid:98) u. A similar analysis shows that the eigenvalues of the matrix R HJ are given by(32) λ ( R HJ ) = ke − ihH ( nx + √ ny ) (cid:16) √ β (cid:16) e ihH ( nx −√ ny ) − e i √ hHny +1 (cid:17) +3 α (cid:16) e i √ hHny (cid:17)(cid:17) h H +6 αk +2 √ βk . In the case of the two triangular meshes, there are two degrees of freedom per generating pattern,corresponding to the elements K j and (cid:102) K j in the generating pattern G j . We write U j = ( u j, , u j, ),where u j, is the average of the solution over the element K j , and u j, is the average of the solutionover (cid:102) K j . The planar wave solution is then given by U j = e i ( n x x j + n y y j ) (cid:98) U , for (cid:98) U = ( (cid:98) u , (cid:98) u ). Weconsider the case of a right-triangular mesh, where the two equal sides of the isosceles right triangleshave length h R = √ √ h E . The method then reads:(33) ∂ t (cid:98) u (cid:98) u = − h R α (cid:98) u − e − ih R n x α (cid:98) u α (cid:98) u + ( β − α ) (cid:98) u − e − ih R n y β (cid:98) u . In the case of the mesh consisting of equilateral triangles, each with side length h E , the methodreads:(34) ∂ t (cid:98) u (cid:98) u = − √ h E (cid:16) √ α + β (cid:17) (cid:98) u + (cid:16) e − ih E n x (cid:16) − √ α + β (cid:17) − e − ih E n y β (cid:17) (cid:98) u (cid:16) − √ α − β (cid:17) (cid:98) u + (cid:16) √ α + β (cid:17) (cid:98) u . omputing the eigenvalues of the corresponding Jacobi matrices R RJ and R EJ , we obtain λ ( R RJ ) = ± ke − ih R ( n x + n y ) √ α (cid:112) β + ( α − β ) e ih R n y h R + 2 αk , (35) λ ( R EJ ) = ± k (cid:0) α + √ β (cid:1) (cid:113) βe ih E n x + (cid:0) √ α − β (cid:1) e ih E n y (cid:0) h E + 6 αk + 2 √ βk (cid:1) (cid:113)(cid:0) √ α + β (cid:1) e ih E ( n x + n y ) . (36)Then, equations (30), (32), (35), and (36) completely determine the speed of convergence forJacobi’s method of each of the four generating patterns considered. In the case of a higher-orderdiscontinuous Galerkin method with basis consisting of piecewise polynomials of degree p > (cid:99) R J j , D , M j , and (cid:98) L j are ( p +1)( p +2)2 × ( p +1)( p +2)2 blocks. In this case, we do not obtain closed-form expressions for theeigenvalues, but rather compute them numerically.We normalize the velocity magnitude and consider β = (cos( θ ) , sin( θ )). On the square mesh, θ can range from 0 to π/
2. On the hexagonal and equilateral triangle meshes, θ ranges from 0 to π/ θ ranges from 0 to π/
4. We consider a fixed spatial resolution h ,and compare the efficiency of the four patterns for three choices of temporal resolution. We firstconsider an “explicit” time step, satisfying the CFL-type condition(37) k exp = h | β | . As one advantage of using an implicit method is that we are not limited by an explicit timesteprestriction of the form (37), we consider three implicit time steps given by k = 3 k exp , k = 2 k ,and k = 4 k . We then maximize over a discrete sample of θ ∈ [0 , π/
4] and over all wavenumbers( n x , n y ), in order to compute maximum eigenvalue for each of the generating patterns. As thenumber of iterations required to converge to a given tolerance scales like the reciprocal of thelogarithm of the spectral radius, we compare the efficiency of the generating patterns by consideringthe ratio log (cid:0) λ max ( R min J ) (cid:1) log ( λ max ( R ∗ J )) , where λ max ( R ∗ J ) is the largest eigenvalue of R ∗ J , for ∗ = H, S, R, E , and λ max ( R min J ) is the smallestamong all λ max ( R ∗ J ). This ratio corresponds to the ratio of iterations required to converge to agiven tolerance when compared with the most efficient among the generating patterns. The resultsobtained for p = 0 , , ,
3, and k = k , k , k for each generating pattern are shown in Table 1 andFigure 2.We remark that for degrees 0, 1, and 2 polynomials, the hexagonal mesh resulted in the small-est eigenvalues for all choices of timestep considered, and the square mesh resulted in the second-smallest eigenvalues. For degree 3 polynomials, the square mesh resulted in the smallest eigenvaluesfor all cases considered. We notice a significant decrease in the expected performance of the hexag-onal elements in the case of p = 3, although we have noticed that the effect observed in practice isnot as significant as the theoretical results would suggest. = 0 p = 1 k k k k k k Squares 1.128939 1.133989 1.136772 1.058098 1.118222 1.130101Right triangles 1.128939 1.133989 1.136772 1.084223 1.132326 1.137313Equilateral triangles 1.207328 1.215467 1.219948 1.137267 1.201638 1.214376 p = 2 p = 3 k k k k k k Right triangles 1.111863 1.126951 1.133634 1.010482 1.005391 1.002733Equilateral triangles 1.177503 1.201918 1.213527 1.074570 1.074570 1.074570
Table 1: Ratio of logarithm of eigenvalues log (cid:0) λ max ( R min J ) (cid:1) / log ( λ max ( R ∗ J )) ranging over angle θ and wavenumber ( n x , n y ), for piecewise polynomials of degree 0, 1, 2, and 3, for varying choices oftime step k . The smallest eigenvalue in each column is highlighted. k k k . . (a) p = 0 k k k . . (b) p = 1 k k k . . (c) p = 2 k k k . . (d) p = 3 Hexagons Squares Right Triangles Equilateral TrianglesFig. 2: Ratios of the logarithm of the largest eigenvalues for each pattern. . Numerical Results.4.1. Advection with variable velocity field. To perform numerical experiments extendingthe analysis of equation (3) beyond the case of a constant velocity β , we consider a variable velocityfield β ( x, y ). In this case, the upwind numerical flux(38) (cid:98) F ( u + , u − , n , x, y ) = (cid:40) u − ( x, y ) if β ( x, y ) · n ≥ u + ( x, y ) if β ( x, y ) · n < β ( x, y ) = (2 y − , − x + 1) on the spatial domain Ω = [0 , × [0 , x , y ) = (0 . , . u ( x, y ) = exp( − x − x ) + ( y − y ) )) . The exact solution is periodic with period π , and is given by the rotation about the center of thedomain,(40) u ( x, y, t ) = exp( − x − . .
15 cos(2 t )) + ( y − . − . t ) sin( t )) )) . . . . . . . . . . . . . Fig. 3: Velocity field β ( x, y ) = (2 y − , − x + 1) We consider meshes of the domaincreated by repeating each of the four generating patterns considered in the previous section. Asbefore, for fixed spatial resolution h , we choose h H , h S , h R , and h E such that the number of degreesof freedom is the same for each mesh. We then solve the advection equation using the backwardEuler time discretization, where the block Jacobi iterative method is used to solve the resulting inear system. The zero vector is used as the starting vector for the block Jacobi solver. We choose h = 0 .
05, and since max ( x,y ) | β ( x, y ) | = √
2, we consider time steps of k = h/ √ k = 2 k , k = 4 k . The number of iterations required for the block Jacobi method to converge to a toleranceof 10 − are given in Table 2. p = 0 p = 1 p = 2 p = 3 k k k k k k k k k k k k
33 57 104 21 41 77
41 77 21 39 75
Squares 35 61 109
42 83
42 83 22 42 81Right triangles 39 68 128 26 51 100 25 51 100 25 51 100Equilateral triangles 37 67 123 25 47 92 25 47 92 24 47 91
Table 2: Iterations required for the block Jacobi iterative method to converge in the case of anon-constant velocity field. The smallest number of iterations in each column is highlighted.The results are similar to those from the analysis performed in Section 3.3. We note that thehexagonal and square meshes resulted in the lowest number of Jacobi iterations for all of the testcases considered. In contrast to the results of Section 3.3, we do not observe a decrease in theperformance of the hexagonal elements for the case of p = 3, and instead the performance is similaramong all choices of p considered. We now consider the effect of polygonal elements onirregular meshes. To this end, we consider a set of generating points distributed evenly on aCartesian grid with mesh size h . Then, each point is perturbed by a random perturbation sampleduniformly from the interval [ − δ, δ ]. We obtain two randomized meshes by constructing the Delaunaytriangulation and Voronoi diagram resulting from this set of generating points. The Delaunay meshconsists entirely of triangular elements, whereas the Voronoi diagram is constructed out of arbitrarypolygonal elements. Examples of the two meshes considered are shown in Figure 4. In contrast tothe regular meshes considered in the previous examples, these two meshes do not consist of the samenumber of elements. The Voronoi diagram consists of about half the number of elements as theDelaunay triangulation. In the test case considered, the randomized polygonal mesh consists of 410polygonal elements, whereas the randomized triangular mesh consists of 759 triangular elements.The governing equations and set-up is the same as in the previous section. We record thenumber of block Jacobi iterations required to converge to a tolerance of 10 − in Table 3. Becausethere is a difference in the number of mesh elements, the resulting linear system will have a differenttotal number of degrees of freedom. This difference will then have an additional effect on the speedof convergence of the block Jacobi method. We note that for polynomials of degree p = 0 , , , k considered, solving the system resulting from the Voronoi diagramrequires fewer block Jacobi iterations than does solving the system resulting from the correspondingDelaunay triangulation. The above analysis focused on the blockJacobi method largely because of the simplicity of the method. In practice, more sophisticatediterative methods are often used [25]. In this section, we consider the solution of the linear system (7)by means of the GMRES method, using both the block Jacobi and the block ILU(0) preconditioners.Since the computational work increases per iteration in GMRES, we choose a restart parameter of20 iterations [29]. We repeat the above test case of the advection equation with variable velocity . . . . . . . . (a) Delaunay triangulation . . . . . . . . (b) Voronoi diagram Fig. 4: Randomized polygonal and triangular meshes corresponding to the same set of generatingpoints. p = 0 p = 1 p = 2 p = 3 k k k k k k k k k k k k
27 32 38 24 33 38 24 32 36 22 31 36
Delaunay triangulation 38 48 52 33 45 48 33 46 50 33 44 48
Table 3: Iterations required for the block Jacobi iterative method to converge in the case of irregular,randomly perturbed meshes. The smallest number of iterations in each column is highlighted.field and record the number of GMRES iterations required to converge to a tolerance of 10 − usingthe block Jacobi preconditioner in Table 4.We now consider the solution of the above problem using the GMRES method with the blockILU(0) preconditioner. Because of the sensitivity of the block ILU(0) factorization to the orderingof the mesh elements, and for the sake of a fair comparison between the generating patterns, weconsider the natural ordering of mesh elements, illustrated in Figure 5. As in the case of the blockJacobi preconditioner, we repeat the test case of the advection equation with variable velocity field.We record the number of GMRES iterations required to converge to the above tolerance using theblock ILU(0) preconditioner in Table 5. In this case, the square mesh resulting in the smallestnumber of iterations in all of the trials. The mesh consisting of right isosceles triangles resulted inthe largest number of iterations in all trials. We further note that the number of GMRES iterationsrequired when using the block Jacobi preconditioner scales similarly to the number of block Jacobiiterations required, as recorded in Table 2. We note that the block ILU(0) preconditioner requiresfewer GMRES iterations to converge, and the number of iterations scales more favorably in k , when = 0 p = 1 p = 2 p = 3 k k k k k k k k k k k k
31 53 92 25 42 80
47 86
49 90
Squares 37 64 116 27 51 101
51 98
52 100Right triangles 40 70 134 33 61 123 31 60 117 29 59 115Equilateral triangles 39 67 124 33 58 113 32 59 113 31 57 111
Table 4: Iterations required for the GMRES iterative method with block Jacobi preconditioner toconverge. The smallest number of iterations in each column is highlighted.compared with the block Jacobi preconditioner.1 2 34 5 67 8 9 (a) Hexagonal mesh (b) Square mesh (c) Right triangular mesh (d) Equilateral triangular mesh
Fig. 5: Illustration of the natural ordering of mesh elements.
The compressible Euler equations of gas dynamics intwo dimensions (see e.g. [14]) are given by(41) u t + ∇ · f ( u ) = 0 , = 0 p = 1 p = 2 p = 3 k k k k k k k k k k k k
10 13 20 11 15 23 10 13 22Squares
Right triangles 13 19 32 10 14 28 10 15 27 11 14 28Equilateral triangles 11 15 27 10 12 22 9 12 22 9 12 22
Table 5: Iterations required for the GMRES iterative method with ILU(0) preconditioner to con-verge. The smallest number of iterations in each column is highlighted.for(42) u = ρρuρvρE , f ( u ) = ρuρu + pρuvρHu , f ( u ) = ρvρuvρv + pρHv , where ρ is the density, v = ( u, v ) is the fluid velocity, p is the pressure, and E is the specific energy.The total enthalpy H is given by(43) H = E + pρ , and the pressure is determined by the equation of state(44) p = ( γ − ρ (cid:18) E − v (cid:19) , where γ = c p /c v is the ratio of specific heat capacities at constant pressure and constant volume.We consider the model problem of an unsteady compressible vortex in a rectangular domain [32].The domain is taken to be a 20 ×
15 rectangle and the vortex is initially centered at ( x , y ) = (5 , θ . The exact solution is given by u = u ∞ (cid:18) cos( θ ) − (cid:15) (( y − y ) − vt )2 πr c exp (cid:18) f ( x, y, t )2 (cid:19)(cid:19) , (45) u = u ∞ (cid:18) sin( θ ) − (cid:15) (( x − x ) − ut )2 πr c exp (cid:18) f ( x, y, t )2 (cid:19)(cid:19) , (46) ρ = ρ ∞ (cid:18) − (cid:15) ( γ − M ∞ π exp(( f ( x, y, t )) (cid:19) γ − , (47) p = p ∞ (cid:18) − (cid:15) ( γ − M ∞ π exp(( f ( x, y, t )) (cid:19) γγ − , (48)where f ( x, y, t ) = (1 − (( x − x ) − ut ) − (( y − y ) − vt ) ) /r c , M ∞ is the Mach number, u ∞ , ρ ∞ , and p ∞ are the free-stream velocity, density, and pressure, respectively. The free-stream velocity isgiven by ( u, v ) = u ∞ (cos( θ ) , sin( θ )). The strength of the vortex is given by (cid:15) , and its size is r c . Wechoose the parameters to be γ = 1 . M ∞ = 0 . u ∞ = 1, θ = arctan(1 / (cid:15) = 0 .
3, and r c = 1 . n the discontinuous Galerkin discretization of the Euler equations we use the Lax-Friedrichsnumerical flux defined by(49) (cid:98) F ( u + , u − , n ) = (cid:0) f ( u − ) · n + f ( u + ) · n + α ( u − − u + ) (cid:1) , where α is the maximum absolute eigenvalue over u − and u + of the matrix B ( u , n ) defined by(50) B ( u , n ) = J f n + J f n , where J f and J f are the Jacobian matrices of the components of the numerical flux function f defined in equation (42).We use the backward Euler time discretization, but remark that (2) results in a nonlinear setof equations, which are solved using Newton’s method. Each iteration of Newton’s method requiressolving a linear equation of the form (7). We set h = 1, and consider three time steps, k = 0 . h , k = 2 k , k = 4 k . We use piecewise polynomials of degrees p = 0 , , ,
3. Each Newton solverequires between 3 to 8 iterations to converge to within a tolerance of 5 × − . The toleranceused for the linear solvers is the same as in the previous test cases. Each iteration of Newton’s method requires the solutionof a linear system of equations. We solve these systems using the block Jacobi method. We computethe total the number of Jacobi iterations required to complete one solve of Newton’s method, andreport the results in Table 6. We note that for each choice of p and time step k , the hexagonal meshrequired the fewest number of block Jacobi iterations. As in the previous numerical experiments,we do not see a decrease in performance for the hexagonal elements in the case of p = 3. The squaremesh resulted in the second-smallest number of iterations for most of the cases considered, whilethe two configurations of triangles resulted in generally similar numbers of iterations. p = 0 p = 1 p = 2 p = 3 k k k k k k k k k k k k
32 49 78 31 50 83 50 90 158 53 97 171
Squares 34 51 89
54 92 54 99 181 55 105 201Right triangles 37 56 97 41 64 112 58 101 189 59 113 217Equilateral triangles 37 57 95 39 62 113 54 99 179 60 114 215
Table 6: Block Jacobi iterations required per Newton solve of the compressible Euler equations.The lowest number of iterations in each column is highlighted.
We now repeat the above test case, using the GMRES methodto solve the resulting linear systems. We consider both the block Jacobi and block ILU(0) precon-ditioners. We then compute the total number of GMRES iterations required to complete one solveof Newton’s method. As in Section 4.1.3, the ordering of the mesh elements has a significant effecton the effectiveness of the block ILU(0) approximate factorization. For this reason, we use thenatural ordering of elements, depicted in Figure 5. We present the results for the block Jacobipreconditioner in Table 7, and for the block ILU(0) preconditioner in Table 8. With the blockJacobi preconditioner, the hexagonal mesh required the smallest number of iterations for all testcases considered, and the square mesh the second-smallest. In the case of the block ILU(0) pre-conditioner, the square mesh required the fewest number of iterations, with the hexagonal mesh sually requiring the second-smallest number of iterations. As we observed in Section 4.1.3, thenumber of iterations required for both the block Jacobi method and GMRES with the block Jacobipreconditioner scales quite poorly with increasing timesteps. The number of GMRES iterationsrequired when using the block ILU(0) preconditioner is significantly better. p = 0 p = 1 p = 2 p = 3 k k k k k k k k k k k k
55 74 106 50 92 126 61 110 153 76 141 195
Squares 62 84 155 52 93 132 67 126 185 78 149 222Right triangles 63 87 162 81 106 184 96 132 242 85 159 299Equilateral triangles 66 90 167 81 108 187 72 133 197 85 161 245
Table 7: GMRES with block Jacobi preconditioner. Iterations required per Newton solve of thecompressible Euler equations. The lowest number of iterations in each column is highlighted. p = 0 p = 1 p = 2 p = 3 k k k k k k k k k k k k
42 21
36 48 29 48 57 29 50 64Squares
24 28
21 33 40 24 41 49 27 48 60
Right triangles 31 40 70 35 40 60 36 48 69 31 49 75Equilateral triangles 28 37 65 37 44 70 33 56 68 38 64 80
Table 8: GMRES with block ILU(0) preconditioner. Iterations required per Newton solve of thecompressible Euler equations. The lowest number of iterations in each column is highlighted.
The following two numerical experiments extend the aboveresults to larger-scale, more realistic flow problems. These problems, in contrast to the precedingtest cases, are characterized by a large number of degrees of freedom, the presence of geometricfeatures and wall boundary conditions, variably-sized mesh elements, and shocks. As in the previoussection, the equations considered here are the compressible Euler equations. For the following twoproblems, we choose the finite element function space to consist of piecewise constant functions(corresponding to p = 0), which results in a finite-volume-type discretization. This choice ofdiscretization allows for the solution of problems with shocks, without the use of slope limiters,artificial viscosity, or other shock-capturing techniques [17]. The Roe numerical flux is used as anapproximate Riemann solver for these problems. For a first test case, we consider the inviscidflow over a circular cylinder at Mach 0.2. The computational domain is defined as Ω = R \ C , where R = [ − , × [ − , C is a disk of radius 1 centered at the point (5 , ∂R , and a no normal flow condition is enforced on ∂C . The freestreamvelocity is taken to be unity in the x -direction, and ρ ∞ = 1. For this test case we use fourunstructured meshes, two consisting entirely of triangles, and two consisting of mixed polygons,generated using the PolyMesher algorithm [31]. All the meshes are created using a gradient-limitedelement size function that determines the initial distribution of seed points according to the rejection
10 0 10 20 30 − − . . . . . . . . . . . . . . Fig. 6: Overview of the coarse mesh with 15,404 elements, with zoom-in showing polygonal elementsnear the surface of the cylinder.method [24], such that the element edge length near the surface of the cylinder is about one-fifththe edge length of elements away from the cylinder. For both the triangular and polygonal meshes,we consider a coarse mesh, with 15,404 elements, and a fine mesh with 62,270 elements. Thus, theaverage area of each element is the same for both the polygonal and triangular meshes. Additionally,the number of degrees of freedom in the solution is the same, allowing for a fair comparison. Thecoarse polygonal mesh, and a zoom-in around the surface of the cylinder are shown in Figure 6.Starting from freestream initial conditions, we integrate the equations until t = 5 × − inorder to obtain a representative solution. Using this solution, we then compute 10 time steps usinga third-order A -stable DIRK method [1]. Each stage of the DIRK method requires the solution ofa nonlinear system of equations, which we solve by means of Newton’s method. In each iterationof Newton’s method, we solve the resulting linear system of the form (7) using both the blockJacobi method and the preconditioned GMRES method. The nonlinear system is solved to withina tolerance of 10 − , and each linear system is solved using a relative tolerance of 10 − . For theGMRES method, we consider two preconditioners: block Jacobi, and block ILU(0). In order tocompare the iterative solver performance differences between meshes, we compute the total numberof solver iterations required to complete all 10 time steps. The results for the GMRES method areshown in Table 9, and for the block Jacobi solver in Table 10.These results demonstrate a consistent trend, corroborating both the numerical results andthe analysis from the previous sections. When using the block Jacobi solver or GMRES with blockJacobi preconditioner, the polygonal mesh results in convergence in between 60–70% of the iterationsrequired for the triangular mesh. The effect is smaller when using the ILU(0) preconditioner, butwe do still observe a modest reduction in the number of iterations required. When using the blockJacobi iterative solver, we observe iteration counts very similar to when using GMRES with blockJacobi as a preconditioner. In these cases, the polygonal mesh requires between 70–80% of theiterations as the all-triangular mesh. able 9: Total GMRES iterations per 10 time steps for inviscid flow over a circular cylinder. (a) Coarse grid with 15,404 elements ILU Jacobi Ratios∆ t Polygonal Triangular Polygonal Triangular ILU Jacobi1 . × −
793 932 2092 3126 0.85 0.672 . × − . × − . (b) Fine mesh with 95,932 elements ILU Jacobi Ratios∆ t Polygonal Triangular Polygonal Triangular ILU Jacobi1 . × − . × − . × − . The next numerical example is designedto investigate the performance of the iterative solvers for steady-state problems, in the presenceof shocks and h -adapted meshes. For this problem, we let the domain be Ω = R \ C , where R = [0 , × [0 , C is a circle of radius one centered at (5 , M = 2 .
0, resulting in the formation of a shock upstream from the cylinder. Inorder to accurately capture the shock, we refine the mesh in its vicinity. As in the previous case,we consider a set of four meshes, two all-triangular, and two polygonal. For both the triangularand polygonal meshes, we consider coarse and fine versions, with 31,162 and 95,932 elements,respectively. The coarse mesh is depicted in Figure 7a, with Mach isolines overlaid to indicate theposition of the shock. Additionally, Mach contours of the steady-state solution are shown in Figure7b. Beginning with freestream initial conditions, the solution rapidly approaches a steady state. We able 10: Total block Jacobi iterations per 10 time steps for inviscid flow over a circular cylinder. (a) Coarse grid with 15,404 elements ∆ t Polygonal Triangular Ratio1 . × − . × − . × − . (b) Fine mesh with 95,932 elements ∆ t Polygonal Triangular Ratio1 . × − . × − . × − . t = 100 in order to obtain an solution which can be used as an initial guessfor the steady-state Newton solve. Then, starting with this solution, we set the time-derivative ofthe solution to zero and solve the resulting nonlinear equations using Newton’s method to find asteady-state solution. The resulting linear system that is required to be solved at each iteration canbe thought of as corresponding to equation (7), where formally we set k = ∞ . The nonlinear systemis solved to within a tolerance of 10 − , and each linear system is solved using a relative toleranceof 10 − . Since the mass matrix in (7) acts to regularize the linear system, the conditioning becomesworse for larger values of k , and the number of iterations required per linear solve grows. Hence,effective preconditioners are particularly important for the solution of such steady-state problems.For these problems, the block Jacobi iterative solver did not converge in fewer than 10,000 iterations,and so we consider only the GMRES method, using block ILU(0) and block Jacobi preconditioners.We present the comparison of iteration counts for this problem in Table 11. On the coarsemeshes, the ILU(0) preconditioner required about 73% as many iterations on the polygonal meshwhen compared with the triangular mesh. This difference is more significant when using the blockJacobi preconditioner, consistent with the results observed in previous section. In this case, thepolygonal mesh requires only slightly more than one third the number of iterations as the all-triangular mesh. On the fine mesh, there are close to half a million degrees of freedom. For aproblem of this scale, we did not observe convergence in less than 10,000 iterations per linearsolve using the block Jacobi preconditioner, and so we only compare performance using the blockILU(0) preconditioner. In this case, the polygonal mesh required about half as many iterations persteady-state solve when compared with the all-triangular mesh.
5. Conclusions.
In this paper we have analyzed the effect of the generating pattern of aregular mesh on the convergence of iterative linear solvers applied to implicit discontinuous Galerkindiscretizations. We considered four generating patters: a hexagon, a square, two right triangles,and two equilateral triangles.A classical von Neumann analysis applied to the constant-velocity advection equation allowedus to compute the eigenvalues of the block Jacobi matrix, and therefore estimate the speed ofconvergence of the block Jacobi method. In more than half of the cases considered, the hexagonalgenerating pattern resulted in the smallest eigenvalues, and in the remaining cases, the square (a) Coarse mesh for supersonic test problem, showingMach isolines for steady-state solution (b) Contours of Mach number for steady state solu-tion Fig. 7: Overview of coarse polygonal mesh with 31,162 elements, showing Mach number contoursfor steady-state solution.Table 11: Total GMRES iterations per steady-state solve for supersonic flow over a cylinder. (a) Coarse grid with 31,162 elements
Polygonal Triangular RatioILU 469 640 0.73Jacobi 2340 6464 0.36 (b) Fine mesh with 95,932 elements
Polygonal Triangular RatioILU 953 1947 0.49Jacobi – – –generating pattern resulted in the smallest eigenvalues.In order to extend these results beyond the case of the constant-velocity advection equation,we performed numerical experiments on the variable-velocity advection equation and compressibleEuler equations. In the case of the advection equation, in all but one case the hexagonal meshresulted in the fastest convergence, and in the remaining case the square mesh resulted in thefastest convergence. In the case of the Euler equations, the hexagonal mesh resulted in the fastest onvergence in all test cases.We additionally considered two irregular meshes resulting from the random perturbation of aset of regularly-spaced generating points. We obtain a triangular mesh by performing the Delaunaytriangulation on these points, and we obtain a polygonal mesh by constructing the Voronoi diagramdual to the Delaunay triangulation. Solving the advection equation on these irregular meshes,we observed that the block Jacobi method converged faster on the polygonal mesh in every testcase. Additionally, we performed numerical experiments examining the performance of the GMRESiterative method when used with the ILU(0) preconditioner. We found that in all of the test cases,the square generating pattern resulted in the fewest number of GMRES iterations, and in all buttwo cases, the hexagonal generating pattern resulted in the second-fewest number of iterations.For a final set of numerical experiments, we performed two inviscid fluid flow simulations onsets of coarse and fine meshes. Each mesh was either all-triangular, or was composed of arbitrarypolygons. We measured iteration counts for both time-dependent and steady-state problems, usingthe block Jacobi method, and GMRES with block ILU(0) and block Jacobi preconditioners. Wefound that the polygonal meshes resulted in faster convergence of the iterative solvers, with a largerdifference being observed for the block Jacobi method and preconditioner. This difference was morepronounced for the steady-state problem, with a quite significant difference observed on the finemesh using GMRES with ILU(0).These results suggest that certain types of polygonal meshes have the advantage of rapid conver-gence of iterative solvers. Future research directions involve the study of accuracy of DG methods onpolygonal and polyhedral meshes, efficient computation of quadrature rules over arbitrary polygonaldomains and the extension of the above results to three spatial dimensions. REFERENCES[1] Roger Alexander. Diagonally implicit Runge-Kutta methods for stiff O.D.E.’s.
SIAM Journal on NumericalAnalysis , 14(6):1006–1021, 1977.[2] Georgios Balafas.
Polyhedral Mesh Generation for CFD-Analysis of Complex Structures . PhD thesis, MasterThesis, Technische Universit¨at M¨unchen, 2014.[3] F. Bassi and S. Rebay.
GMRES Discontinuous Galerkin Solution of the Compressible Navier-Stokes Equations ,pages 197–208. Springer Berlin Heidelberg, Berlin, Heidelberg, 2000.[4] Michele Benzi, Wayne Joubert, and Gabriel Mateescu. Numerical experiments with parallel orderings for ILUpreconditioners.
Electronic Transactions on Numerical Analysis , 8:88–114, 1999.[5] Martin Berggren. A vertex-centered, dual discontinuous Galerkin method.
J. Comput. Appl. Math. , 192(1):175–181, 2006.[6] Bernardo Cockburn and Chi-Wang Shu. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems.
J. Sci. Comput. , 16(3):173–261, 2001.[7] E. Cuthill and J. McKee. Reducing the bandwidth of sparse symmetric matrices. In
ACM Proceedings of the1969 24th National Conference , ACM ’69, pages 157–172, New York, NY, USA, 1969. ACM.[8] Laslo T. Diosady and David L. Darmofal. Preconditioning methods for discontinuous Galerkin solutions of theNavier-Stokes equations.
J. Comput. Phys. , 228(11):3917–3935, June 2009.[9] Boris Diskin and James L Thomas. Comparison of node-centered and cell-centered unstructured finite-volumediscretizations: inviscid fluxes.
AIAA journal , 49(4):836–854, 2011.[10] Boris Diskin, James L Thomas, Eric J Nielsen, Hiroaki Nishikawa, and Jeffery A White. Comparison ofnode-centered and cell-centered unstructured finite-volume discretizations: viscous fluxes.
AIAA journal ,48(7):1326–1338, 2010.[11] Iain S. Duff and G´erard A. Meurant. The effect of ordering on preconditioned conjugate gradients.
BITNumerical Mathematics , 29(4):635–657, 1989.[12] Rao V. Garimella, Jibum Kim, and Markus Berndt. Polyhedral mesh generation and optimization for non-manifold domains. In
Proceedings of the 22nd International Meshing Roundtable , pages 313–330. Springer,2014. 2213] Alan George. Nested dissection of a regular finite element mesh.
SIAM Journal on Numerical Analysis ,10(2):345–363, 1973.[14] Ralf Hartmann. Discontinuous Galerkin methods for compressible flows: higher order accuracy, error esti-mation and adaptivity. In H. Deconinck and M. Ricchiuto, editors,
VKI LS 2006-01: CFD-Higher OrderDiscretization Methods, Nov. 14-18, 2005 . Von Karman Institute for Fluid Dynamics, Rhode Saint Gen`ese,Belgium, 2005.[15] J. S. Hesthaven and T. Warburton.
Nodal discontinuous Galerkin methods , volume 54 of
Texts in AppliedMathematics . Springer, New York, 2008.[16] Ethan J. Kubatko, Clint Dawson, and Joannes J. Westerink. Time step restrictions for Runge-Kutta discon-tinuous Galerkin methods on triangular grids.
J. Comput. Phys. , 227(23):9697–9710, December 2008.[17] Randall J. LeVeque.
Finite volume methods for hyperbolic problems , volume 31. Cambridge university press,2002.[18] Hong Luo, Joseph D. Baum, and Rainald L¨ohner. A discontinuous Galerkin method based on a Taylor basisfor the compressible flows on arbitrary grids.
J. Comput. Phys. , 227(20):8875–8893, 2008.[19] Gianmarco Manzini, Alessandro Russo, and N. Sukumar. New perspectives on polygonal and polyhedral finiteelement methods.
Math. Models Methods Appl. Sci. , 24(8):1665–1699, 2014.[20] Harry M. Markowitz. The elimination form of the inverse and its application to linear programming.
Manage-ment Science , 3(3):255–269, 1957.[21] Wayne Oaks and Stefano Paoletti. Polyhedral mesh generation. In
Proceedings of the 9th International MeshingRoundtable , pages 57–67, 2000.[22] J. Peraire, M. Vahdati, K. Morgan, and O. C. Zienkiewicz. Adaptive remeshing for compressible flow compu-tations.
J. Comput. Phys. , 72(2):449–466, 1987.[23] M. Peric. Flow simulation using control volumes of arbitrary polyhedral shape.
ERCOFTAC bulletin , 62:25–29,2004.[24] Per-Olof Persson.
Mesh generation for implicit geometries . PhD thesis, Massachusetts Institute of Technology,2005.[25] Per-Olof Persson. Scalable parallel Newton-Krylov solvers for discontinuous Galerkin discretizations. In . Amer-ican Institute of Aeronautics and Astronautics, January 2009.[26] Per-Olof Persson and Jaume Peraire. Newton-GMRES preconditioning for discontinuous Galerkin discretiza-tions of the Navier-Stokes equations.
SIAM Journal on Scientific Computing , 30(6):2709–2733, 2008.[27] W. H. Reed and T. R. Hill. Triangular mesh methods for the neutron transport equation. Technical ReportTechnical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973.[28] J. Ruppert. A Delaunay refinement algorithm for quality 2-dimensional mesh generation.
J. Algorithms ,18(3):548–585, 1995.[29] Yousef Saad.
Iterative methods for sparse linear systems . Society for Industrial and Applied Mathematics,Philadelphia, PA, second edition, 2003.[30] Jonathan R. Shewchuk. Delaunay refinement algorithms for triangular mesh generation.
Comput. Geom. ,22(1-3):21–74, 2002.[31] Cameron Talischi, Glaucio H. Paulino, Anderson Pereira, and Ivan F. M. Menezes. Polymesher: a general-purpose mesh generator for polygonal elements written in matlab.
Structural and Multidisciplinary Opti-mization , 45(3):309–328, 2012.[32] Z.J. Wang, Krzysztof Fidkowski, R´emi Abgrall, Francesco Bassi, Doru Caraeni, Andrew Cary, Herman Decon-inck, Ralf Hartmann, Koen Hillewaert, H.T. Huynh, Norbert Kroll, Georg May, Per-Olof Persson, Bramvan Leer, and Miguel Visbal. High-order CFD methods: current status and perspective.