The interior penalty discontinuous Galerkin method for elastic wave propagation: grid dispersion

it’s dispersion particular to the effect that different basis functions have on the dispersion. consider different types of functions that naturally yield a diagonal mass matrix. is to seismology because a diagonal mass matrix is tantamount to an explicit and efﬁcient time marching scheme. the Legendre basis functions that are traditionally used in the introduce numerical dispersion and anisotropy. that using Lagrange basis functions along with the Gauss nodes has attractive advantages for numerical wave propagation.

determines the largest time step for the time-marching scheme such that the numerical solution remains bounded. On the other hand, the grid dispersion criterion determines the largest sampling ratio for the spatial discretization, conveniently denoted as nodes per wavelength, such that the numerical solution has an acceptable accuracy.
Although the FEM has been applied to the wave propagation problem since the late 1960s and early 1970s (Chopra et al. 1969;Lysmer & Drake 1972), the grid dispersion and stability analysis was not available until the early 1980s, and only for the first degree FEM schemes (Mullen & Belytschko 1982;Marfurt 1984). A general framework to analyse higher-degree methods was developed in Cohen (2002). His approach overcomes the difficulties due to irregular grid spacing. He analysed the spectral element method (SEM) for the acoustic case and showed results up to third degree basis functions. The hp-version of the FEM for acoustic wave propagation was analysed in Ainsworth (2004a). In his analysis, he obtained closed form relations for the grid dispersion using Padé approximants. We have recently extended the grid dispersion and stability analysis to the high-degree FEM for the elastic case in De Basabe & Sen (2007), where we considered the continuous and finite-difference cases for the time derivative.
The grid dispersion of the DGM for wave propagation has been analysed in Hu et al. (1999), Stanescu et al. (2000), Ainsworth (2004b) and Ainsworth et al. (2006). In Hu et al. (1999), the dispersion and dissipation errors of the discretization of the scalar advection equation and the acoustic wave equation in one space dimension are considered using the flux formulation DGM, Legendre basis functions, and triangular and quadrilateral elements. Stanescu et al. (2000) considered the flux formulation of the DGM applied to the scalar advection equation and linearized Euler equation in one spatial dimension. They did not use plane wave analysis and considered the dispersion due to the boundary conditions. Their formulation depends on the particular discretization of the domain. Ainsworth (2004b) studied the linearized advection equation in multiple space dimensions using a flux formulation of the DGM and tensor product basis functions. Finally, Ainsworth et al. (2006) considered the acoustic wave equation and the IP-DGM as well as the flux formulation. Their grid dispersion results include up to third degree polynomial basis functions and conjecture on the extension to higher degree.
It should be noted that, to the best of our knowledge, there are no available results for the grid dispersion properties of any of the DGM applied to the elastic wave equation. Furthermore, it cannot be assumed that the results available for the acoustic case can be applied to the elastic case, since that is not the case in the finite differences method or in the continuous FEM, except for low-order time-stepping schemes (De Basabe & Sen 2007).
The goal of this paper is to present the grid dispersion analysis of the IP-DGM as applied to the elastic wave equation in 2-D. For the analysis, we will consider square elements with tensor-product basis functions. The approach presented here can be extended to 3-D and to other types of elements, as long as they are periodic. We will consider three types of basis functions that naturally yield a diagonal mass matrix and describe the effect that each of them has on the dispersive properties of the numerical scheme. The paper is organized as follows. In the next section, we briefly introduce the IP-DGM. In Section 3, we derive the grid dispersion relations; In Section 4, we describe the results, and in Section 5 we summarize the conclusions.

Weak formulations of the elastic wave equation
Here we will introduce the weak (or variational) formulation of the elastic wave equation used in the IP-DGM. In order to do this, we first introduce the strong form of the elastic wave equation in an isotropic, elastic and unbounded domain, which is given by with u(x, t) → 0 for all t as |x| → ∞ and u( is the density, and λ(x) and μ(x) are the Lamé parameters.
Let h be a finite element partition of the domain, and let h be the set of all the faces between the elements in h . It can be easily shown (see Appendix A) that eq. (1) is equivalent to the following interior-penalty weak formulation: τ γ is the traction vector, given in the isotropic case by and n γ is a unit vector normal to γ . In the above equation we use Einstein summation convention and the shorthand notation u i, j = ∂u i /∂ x j . The double dot product in eq. (4) is given by ∇u : ∇v = u i, j v i, j . In eq. (5), we are using {·} to denote the average of the function and [·] to denote the function jump. This is the standard notation in the DGM literature (Wheeler 1978), and it is defined as follows. Let γ be the edge between the elements E 1 and E 2 , then the jump and average of a vector function u on γ are given by The parameter R in eq. (5) is the penalty, and S is a parameter that takes the values 0, 1 or −1 depending on the particular IP-DGM: S = −1 for SIPG (Darlow 1980), S = 1 for NIPG (Rivière et al. 1999 and S = 0 for incomplete interior penalty Galerkin (IIPG, Dawson et al. 2004). For SIPG, J γ is symmetric with respect to u and v, that is, J γ (u, v; S, R) = J γ (v, u; S, R), which together with the symmetry of B E implies that the stiffness matrix is symmetric. This formulation applied to the acoustic wave equation has been found to be optimally convergent in the energy and L 2 norms with respect to the sampling ratio (Grote et al. 2006). On the other hand, for NIPG, J γ satisfies the antisymmetric condition J γ (u, v ; S, 0) = −J γ (v, u ; S, 0). This formulation was studied in Rivière & Wheeler (2003) for acoustic and elastic wave propagation and it was found to be optimally convergent in the energy norm but suboptimal in the L 2 norm. Finally, for IIPG, J γ has a simplified form but does not yield a symmetric stiffness matrix as in the SIPG formulation. This formulation has not been used for wave propagation to the best of our knowledge, but has been shown to offer advantages for example in flow problems (Dawson et al. 2004).
To solve for the displacement in eq. (2), we introduce a subspace of V h by defining a finite number of basis functions φ E i , i = 1, . . . , m in element E for all E ∈ h . The basis functions φ E i will be discussed in the next subsection. Upon substituting the test function and the displacement by linear combinations of the basis functions we obtain a system of ordinary differential equations that can be solved for the displacement.

Basis functions
A general description of the basis functions can be found in classical finite element literature, for example, Hughes (2000). The nodal and modal basis functions for higher-degree methods are described in detail in Karniadakis & Sherwin (2005). The above references deal with basis functions for the FEM, and therefore, define the basis functions to be continuous across the entire domain. A description of the basis functions used in the DGM can be found in Li (2006). An important difference between FEM and DGM is that in the DGM the basis functions are not required to be continuous over the entire domain but only inside the elements. In general, all the basis functions used in the FEM can be used in the DGM with simplifications, since in the DGM the basis functions can be defined locally on each element. This important feature of the DGM implies that the mass matrix is always block-diagonal, which translates into an efficient time marching algorithm. Furthermore, the basis functions can be chosen such that the mass matrix is exactly diagonal. This is a desirable property because it is necessary to invert the mass matrix for the time marching algorithm. In the following we will restrict our attention to tensor product basis functions in quadrilateral elements that yield a diagonal mass matrix. We will use κ to denote the polynomial degree of the basis functions in one side of the element and m = (κ + 1) 2 to denote the total number of basis functions in the elements.
The first approach to obtain a diagonal mass matrix is to use tensor products of the Legendre polynomials as the basis functions. The Legendre polynomials up to degree 3 are shown in Fig. 1(a). These are called modal basis functions (Karniadakis & Sherwin 2005), and it is the traditional approach that has been used in DGM when a diagonal mass matrix is sought Li 2006). They are orthogonal under the L 2 inner product and have simple recursion formulas for the higher degree polynomials and their derivatives. It should be noted that this approach yields a diagonal mass matrix only if the media pa- rameters are constant inside each element, and that the condition number of the mass matrix is greater than 4κ 2 in the 2-D case (see Appendix B). A second approach is to use nodal basis functions (Karniadakis & Sherwin 2005), making the nodes inside the element match the quadrature points. This idea has been exploited in the SEM, where the Gauss-Lobatto-Legendre (GLL) points and quadrature rules are used to impose the continuity of the basis functions at the element boundaries (Cohen et al. 1994;Komatitsch & Vilotte 1998). The third degree Lagrange polynomials using the GLL nodes are shown in Fig. 1(b). This approach leads to a diagonal mass matrix independently of how the media parameters change inside the elements, but the mass matrix integrals are not computed exactly even if the media parameters are piecewise constant.
The third approach that we will consider is closely related to the second one, but using the Gauss nodes and quadrature rules instead of the GLL nodes and quadrature rules. The Lagrange polynomials using these nodes are shown in Fig. 1(c). The main difference between the Gauss and the GLL nodes is that the endpoints of the interval are always included in the GLL nodes. We can use the Gauss nodes in DGM to define the basis functions because, unlike in FEM, the basis functions do not have to be continuous across the elements. In this approach, as well as in the second one, the mass matrix is always diagonal but, unlike in the second approach, the integration will be exact for piecewise constant and piecewise linearly varying media parameters because we gain 2 degrees of accuracy.

P L A N E WAV E A N A LY S I S
In this section, we apply the plane wave analysis to the IP-DGM to investigate it's grid dispersion properties. For the analysis, we assume that the media is isotropic, unbounded, homogeneous and source free; as mentioned before, these assumptions are always made whenever the plane wave analysis is applied, see for example, Alford et al. (1974), Mullen & Belytschko (1982), Marfurt (1984), Hu et al. (1999), Cohen (2002), Ainsworth (2004a,b), Ainsworth et al. (2006) andDe Basabe &Sen (2007). We emphasize that it is not expected that practical applications satisfy these assumptions, they are for analysis purposes only, to understand the basic accuracy of the numerical scheme. Nevertheless, the outcome of the analysis will be useful to set the simulation parameters in realistic applications. Furthermore, we will assume that the finite element mesh is periodic, and that the elements are square with sides parallel to the coordinate axis and with tensor product basis functions; these are common assumptions when a FEM is analysed, see for example, Marfurt (1984), Cohen (2002), Ainsworth (2004a,b), Ainsworth et al. (2006) andDe Basabe &Sen (2007).
Introducing the assumptions into the weak formulation of the elastic wave equation, eq. (2), we obtain where α = (λ + 2μ)/ρ is the P-wave velocity, β = √ μ/ρ is the S-wave velocity and r = α/β is the ratio of the P-to S-wave velocity.
In order to discretize eq. (8) using the IP-DGM, we introduce the following approximation to the displacement vector using the basis functions: where ξ E i and η E i are the coefficients of the x and z components of displacement in element E. We also write the test function using the basis functions as follows where a E i and b E i are arbitrary coefficients. Without loss of generality we can set a E i = 1 for E = E 0 and i = j, a E i = 0 otherwise, where E 0 and j are arbitrary, and b E i = 0, to obtain Setting now a E i = 0 and b E i = 1 for E = E 0 and i = j, b E i = 0 otherwise, and substituting in eq. (8) yields Note in the last term of the left-hand side of eqs (13) and (14) that there are only four elements of h for which J γ is non-zero, we call these γ T , γ B , γ L and γ R (see Fig. 2). Using the linearity of B E and J γ with respect to the first argument and computing the integrals on the master elementÊ, we can write eqs (13) and (14) as and where and whereγ f , f ∈ {T, B, L , R} are the edges of the master elementÊ. If we assume that the displacement is a plane wave, then and where k is the wavenumber, x i contains the ith node coordinates and A i and B i are arbitrary constants. The plane wave assumption implies that and and similar expressions for η E f i , f ∈ {T, B, L , R}. Substituting these in eqs (15) and (16) we obtain the following generalized eigenvalue problem of order 2m: where = h 2 ω 2 h /β 2 , ω h is the angular frequency at which the wave travels in the grid and K ν i j , ν = 1, . . . , 4, are the so-called dynamic stiffness matrices, given by (34) Note that we have not assumed any particular kind of basis functions or grid nodes for the eigenvalue problem in eqs (32) and (33). The basis functions can be any of the ones described in Section 2 or other. The number of eigenvalues will usually exceed the number of physical modes, therefore, we need to identify which eigenvalues correspond to the P and S waves. We can easily do this by calculating the velocities corresponding to each eigenvalue and comparing to the known P-and S-wave velocities (Sármány et al. 2007). Let us denote as P and S the eigenvalues corresponding to the P and S waves.
The grid dispersion of the P and S waves is given by the ratio between the velocity at which the wave travels in the grid and the physical velocity. From the definition of the eigenvalues we have that the angular frequency of the S wave in the grid is given by ω h = (β/ h) √ S , therefore, the velocity at which the S wave travels in the grid is given by where δ = h/(κL) is the sampling ratio and L is the wavelength. It is convenient to measure the grid dispersion as the relative error in the velocity, given by and similarly for the error in the P-wave velocity In the above error expressions, the sign of the error will indicate if the numerical scheme causes the waves to be delayed or hastened. The grid-dispersion error will depend on the sampling ratio δ, the P-to S-wave velocity ratio r, the wavenumber k and the degree of the basis functions κ, as well as on the IP-DGM parameters S and R.

R E S U LT S
We now proceed to discuss the accuracy of the IP-DGM from three different perspectives. We will consider (i) the convergence of the methods with respect to the sampling ratio, (ii) the convergence with respect to the polynomial degree of the basis functions and (iii) the anisotropy introduced by the grid dispersion. We finally compare the methods with the SEM, which is a method that has become very popular for seismic wave propagation. The grid dispersion of the SEM will be computed using the approach of De Basabe & Sen (2007). We will focus on the grid dispersion of the S wave since it is always larger than the dispersion of the P wave and thus it is of more concern. The convergence with respect to the sampling ratio of the IP-DGM of degrees 1-4 using the GLL basis functions, r = 10 (Poisson's ratio of 0.495) and an incidence angle of θ = 45 • is shown in Fig. 3. This type of analysis is very common in the finiteelement literature. Often a polynomial relation between the absolute value of the error and the size of the elements is assumed of the form |e S | = O(h q ), where q is the convergence rate. Clearly the convergence rate is related to the slope of the convergence curves displayed in Fig. 3. As a visual aid, line segments are displayed to indicate the slopes of the convergence curves of different degrees. It is clear from Fig. 3(a) that q ≈ 2κ for the SIPG method using the GLL basis functions. The same convergence rates are achieved for the Gauss basis functions, but that is not the case for the NIPG and IIPG methods and for the Legendre basis functions, in which cases slower convergence rates are observed. Note that the convergence rates for NIPG and IIPG are of order κ + 1 for odd degree and κ for even degree (Sun & Wheeler 2005). The superconvergence of the SIPG formulation and the different convergence rates of even and odd degrees of the IIPG and NIPG formulations are remarkable results that demand further research.  The convergence of the methods with respect to the degree of the basis functions is shown in Fig. 4. For this figure we have used the following parameters: δ = 0.1 (10 gridpoints per wavelength), θ = 45 • , and r = 10. A consistent feature in this figure is that the convergence rate is slower when the Legendre basis functions are used. We conjecture that this is because the condition number of the mass matrix increases rapidly when the degree of the basis functions is increased. On the other hand, the GLL and Gauss basis  functions have faster and similar convergence rates. The convergence rates of the three methods using the Gauss basis are compared in Fig. 4(d). For the SIPG method a maximum accuracy of approximately 10 significant digits is achieved for κ = 4, and after this point no further improvement is observed. For the NIPG and IIPG methods the maximum accuracy is achieved at approximately κ = 5. The anisotropy introduced by the numerical schemes is displayed in Figs 5 and 6. For visualization purposes we have magnified the dispersion by using r = 10, the Legendre basis functions and calculating the dispersion after a propagation of 200 wavelengths. Furthermore, we have used δ = 0.2 (5 gridpoints per wavelength). In this figure we consider κ = 2, 3 and 4 only. The grid dispersion for κ = 1 is much larger, and therefore, this polynomial degree is unsuitable for practical applications (see Fig. 7). On the contrary, it is uninteresting to display the grid dispersion for κ > 4 because it is very small and isotropic for all practical purposes. The dispersion is quite small and isotropic for κ > 2 for all the methods, as can be seen in Figs 5(a), (b) and (c). Note that for κ > 2 the waves are slightly delayed in all the methods and all incidence angles. To emphasize that using the GLL and Gauss basis functions the method exhibits less anisotropy, we show in Fig. 6 the anisotropy curves for the SIPG method using these basis functions. It can be observed in this figure that the grid dispersion is very small and practically isotropic for κ > 2. Similar results are observed for NIPG and IIPG using nodal basis functions of degree greater than 2.
We now proceed to compare the grid dispersion of the IP-DGM with that of the continuous case. The grid-dispersion curves for the first degree methods are shown in Figs 7 and 8. The grid dispersion of the first degree methods is identical for SIPG, NIPG and IIPG, the only difference is whether the mass matrix is consistent or lumped. We observe by comparing Figs 7(a) with 8(a) and 7(b) with 8(b) that the dispersive behaviour of the first degree methods using the Legendre or Gauss basis is similar to that of the FEM with consistent mass matrix, and using the GLL basis yields a grid dispersion similar to that of the SEM. In fact, Figs 8(a) and (b) are the limiting case of using an infinite penalty in Figs 7(a) and (b). A consistent feature in these figures is that large and anisotropic errors are introduced by the grid dispersion in all the first order methods. Note for example that at δ = 0.1 (10 gridpoints per wavelength) the error is up to 50 per cent using the Gauss or Legendre basis and up to 100 per cent using the GLL basis.
We also compare the convergence with respect to the sampling ratio, the convergence with respect to the polynomial degree of the basis functions and the anisotropy introduced by the grid dispersion for the SEM in Fig. 9. A comparison of Fig. 9(a) with Fig. 3(a) and of Fig. 9(b) with Fig. 4(a) reveals that the convergence rates of the SEM method are similar to that of the SIPG method using the GLL or Gauss basis functions. Also, comparing Fig. 9(c) with Fig. 6(a) we observe that the anisotropic distribution of the grid dispersion of the SEM is very similar to that of the SIPG method with GLL basis functions but larger than that of using the Gauss basis functions. For all the above tests we have used a penalty of R = 1000. Our numerical experiments indicate that the grid dispersion does not change significantly as a function of the penalty whenever it is large enough, but numerical anisotropy is introduced in some combinations of methods and basis functions for R < 1000.

C O N C L U S I O N S
We have developed a method to calculate the grid dispersion of the IP-DGM that is general enough to accommodate nodal or modal basis functions and non-equispaced nodes. The limitations of this approach are that it assumes regular quadrilateral elements and it does not take into account the boundary conditions.
Based on our results we conclude that the SIPG method has attractive advantages over the other formulations. Namely, it allows for lower polynomial degree and sampling ratio to be used to get high accuracy. The results indicate that high accuracy and isotropy can be obtained using the SIPG method of degree 4 or greater with 4 gridpoints per wavelength (δ = 0.25). We have also compared the effect that different basis functions have on the accuracy of the methods. Based on an accuracy criterion, nodal basis functions using the GLL or Gauss nodes have attractive advantages over modal basis functions using the Legendre polynomials. However, we remark that the Legendre basis functions have the advantage of being hierarchical; this may be an important characteristic in some applications to achieve a better performance.
In order to put the results in perspective, we have compared the accuracy of the IP-DGM with that of the SEM. As expected, we have found that the SIPG method with the GLL basis functions performs with practically the same accuracy as the SEM. In practical applications where the media parameters change inside the elements, higher accuracy than the SEM can be expected if the Gauss basis functions are used, since in those cases the under integration of the SEM will deteriorate the accuracy whereas using the Gauss basis functions and quadrature rules will give better approximations. It should be noted that the IP-DGM has the further advantage over the SEM that it can handle non-conforming finite-element meshes (e.g. Sun & Wheeler 2005;.
We have focused on the interior penalty formulation of the DGM, therefore the conclusions are not applicable to the flux formulations, as for example the ADER-DG method . However, the basis functions considered here can be used in a flux formulation, therefore we expect that a grid dispersion study of this type of methods and comparing the basis functions considered here should arrive to the same conclusions. Further research is needed to ascertain this hypothesis.
We have analysed the grid dispersion in the semi-discrete case only. The discretization of the time derivative is an important source of dispersion also, therefore a high-order discretization should be used to preserve the accuracy of the method (see Mercerat et al. 2006). Other sources of grid dispersion that have not been studied here include heterogeneities in the media parameters and irregularities in the finite-element mesh (see Cohen 2002).

A C K N O W L E D G M E N T S
This work was done while the first author was supported by a PhD fellowship from the Mexican National Council for Science and Technology (CONACYT) and a PhD summer support by the EDGER forum and UTIG fellowship.
where n is the quadrature order, and x k and w k are the quadrature nodes and weights. Recall that the basis functions are defined using the Lagrange polynomials and the quadrature nodes, and that for the GLL basis the nodes and weights are those of the GLL quadrature rules, whereas for the Gauss basis those of the Gauss quadrature rules. The condition number of the mass matrix for the nodal basis is therefore given by cond (M ii ) = max k (w k )/ min k (w k ).
It should be noted that the condition number does not create any problems in the implementation of these basis functions because the inversion of a diagonal matrix is a stable procedure. The only concern is the accuracy.