Some comparisons and analyses of time or space discontinuous Galerkin methods applied to elastic wave propagation in anisotropic and heterogeneous media

This research work presents some comparisons and analyses of the time discontinuous space–time Galerkin method and the space discontinuous Galerkin method applied to elastic wave propagation in anisotropic and heterogeneous media. Mechanism of both methods to ensure their stability using time or space discontinuities of unknown fields is analyzed and compared. The most general case of anisotropic and heterogeneous media with physical interfaces of discontinuous material properties is considered, especially for the space discontinuous Galerkin method. A new stability result is proved. Numerical applications to different elastic media, more particularly polycrystalline materials containing a large number of physical interfaces, are also presented to confirm theoretical analyses.


Introduction
Nowadays, as a classical problem, the numerical modeling of elastic wave propagation can be done reliably and efficiently in a large number of cases. However, for heterogeneous media and when the domain of propagation is large compared to the involved wavelengths and to the characteristic length of heterogeneities, accurate and efficient simulations of wave propagation phenomena still remain a challenging task. The case where the geometries of interior physical interfaces of heterogeneities are complicated, for instance the polycrystalline materials, and the 3D case are even more problematic. Therefore, the research field for the development of effective solvers is still active and it is also important to compare and assess the existing methods that have already proved their effectiveness in simpler cases.
The purpose of the paper is to present some comparisons and analyses between a space discontinuous Galerkin (dG) method and a time discontinuous space-time Galerkin method, hereafter named the time dG method. Both methods use discontinuities of unknown fields to ensure the stability but in different ways leading to different consequences.
The time dG method considered herein is a displacement and velocity two-fields method, which is a specific case of a general time discontinuous space-time finite element method for the second order hyperbolic elastodynamic equations developed by Hughes and Hulbert [1][2][3]. In the weak formulation of the general method, additional stabilizing terms of least-squares type are introduced so it can be viewed as a Petrov-Galerkin method. However, these least-squares terms are not necessary for the stability of the method and can be omitted. Without these terms, we get a time discontinuous Galerkin method. In this case, the only appropriate definition of upwind fluxes (due to the causality), together with the use of the elastic energy inner product to enforce the displacement-velocity compatibility, is sufficient to make the time dG method stable (see "Displacement-velocity two fields time dG method" section). By choosing specific spacetime elements, the time dG method can be finally written as a time-stepping scheme and results in an implicit solver. As one important advantage, the method provides an appropriate framework to develop adaptive computing as it remains unconditionally stable even if the space discretization changes over time [4][5][6][7][8]. From one time step to another, the energies that are dissipated in the jumps in time of both displacement and velocity fields guarantee the unconditional stability.
The space dG method 1 is based on the use of spatially element-wise discontinuous finite element basis functions. It has been widely used for computational fluid dynamics, and its application to solid elastodynamics is more recent. We refer to [9][10][11][12][13][14][15][16][17][18] for a non-exhaustive review of related research works. To apply the space dG method, the second-order elastic wave equations are transformed to a first-order hyperbolic system, both velocity-stress [10,11] and velocity-strain [15] formulations can be used. One important advantage of the space dG method is that the use of discontinuous finite element basis functions allows element-wise solving and makes straightforward the parallel data structures and computing.
For the space dG method, developing appropriate numerical fluxes on element interfaces is the key point for its success. For this purpose, exact solving of the Riemann problem defined on element interfaces is usually recommended. Herein, the upwind numerical fluxes that are exact solutions of the Riemann problem are called exact upwind numerical fluxes and those are only approximate solutions are called approximate upwind numerical fluxes. When the space dG method is applied to elastic media, due to the involvement of the fourth order elastic tensor, the exact upwind numerical fluxes can be easily defined only in the case of continuous material properties [11,15]. In [15], the exact upwind numerical fluxes were also developed for 2D isotropic elastic media with discontinuous material properties for the velocity-strain formulation.
Recently, in the most general case of multidimensional anisotropic elastic media with discontinuous material properties, a unified and wave oriented variational framework has been proposed by the author [17], which allows a systematic development of the exact upwind numerical fluxes for both velocity-stress and velocity-strain formulations [18]. Owing to the compact and intrinsic forms expressed in terms of tensors, the derivation of upwind numerical fluxes can be done in a well structured way, which allows better understanding of the physical meaning of the developed upwind numerical fluxes. The consistency of the exact numerical fluxes can be easily proved [15,18]. However, the stability analysis of the space dG method is much more difficult. Wilcox et al. has proved the stability of the velocity-strain dG method in 2D isotropic elastic media including physical interfaces. Their result shows that the energies dissipated in the jumps in space of the velocity and strains fields guarantee the stability for the velocity-strain dG formulation that is not yet discretized in time [15].
In the present paper, the velocity-stress space dG formulation is considered in the general case of anisotropic elastic media with physical interfaces and its stability is analyzed. The result allows explaining the following phenomena previously observed in [17]: when the degree of discontinuity at physical interfaces is high, the use of approximate numerical fluxes leads to instability problems, which cannot be resolved by simply reducing the time step used by the time discretization scheme.
Several numerical examples are presented to illustrate the analytical results. More particularly, polycrystalline materials that present an interesting case of heterogeneous media with a large number of physical interfaces are considered. For this, ultrasonic wave propagation in a single-phase and untextured polycristal composed of elliptic grains is simulated.

Displacement-velocity two fields time dG method
We consider the wave propagation in an elastic medium ⊂ R d of space dimension d (d = 1, 2, 3) and in a time interval [0, T ]. is submitted to dynamic body and surface forces f and g. g is applied on a part of 's boundary ∂ N ⊂ ∂ .
Within the framework of the displacement-velocity two fields time dG method, the elastic wave propagation problem is governed by the following first-order in time equations with as primary unknowns the displacement field u(x, t) and the velocity field v(x, t): where ρ denotes the density and σ the second-order Cauchy stress tensor. σ is linked to the second-order infinitesimal strain tensor ε and the displacement u by the Hooke's law and the definition of ε under the hypothesis of small deformations: In (2), C is the fourth-order elasticity Hooke tensor, (·) T denotes the transpose operator and ":" the usual double dot product between two tensors defined as (C : ε) ij = k,l C ijkl ε kl ≡ C ijkl ε kl . Herein, the Einstein summation convention is systematically used and all the vectors and tensors are denoted using bold letters. The partial derivative operator with respect to time is denoted ∂ t .
Otherwise, it is useful to recall the usual space gradient and divergence operators defined using an orthonormal basis (e i ) i=1,...,d : with "⊗" the usual tensor product between two vectors: (a⊗b) ij = a i b j and "·" the usual dot product between a tensor and a vector: (A · a) i = A ij a j . In the following, the symmetrized tensor product "⊗ s " will also be used and it is defined as: (a ⊗ s b) = 1 2 (a ⊗ b + b ⊗ a). Among the two equations of the governing system (1), the first one (1a) expresses the elastodynamic equilibrium and the second one (1b) the compatibility condition between the displacement and velocity fields, which are treated as independent unknowns. However, (1b) is less classical due to the use of the space derivative operator Div x (σ(·)), which allows, as will be indicated later, proving the stability of the time dG method in terms of kinetic and elastic energies. It is worth noticing that, due to the special form of the equation (1b), the extension of the time dG method to nonlinear problems is not straightforward, since the meaning of the term Div x (σ(∂ t u − v)) in (1b) becomes unclear. For example, in the case of an elastoplastic problem, if only the elastic part of the stress operator is taken into account to have (1b) well-defined, then the stability analysis of the obtained solver would be difficult.
For the time-dependent problem (1), the following initial conditions are necessary: with u 0 and v 0 respectively the initial displacement and velocity fields. Finally, to complete the definition of the elastic wave propagation problem, the following boundary conditions are prescribed: The main idea of the displacement-velocity two fields time dG method is to establish a space-time variational formulation of the strong formulation (1), which is furthermore already discretized in time.
To do this, the whole space-time domain S = × ]0, T [ is subdivided into N spacetime slabs: {S n = × T n } n=1,···,N , with T n = ]t n−1 , t n [. Between two successive spacetime slabs, the unknown fields u and v can be discontinuous. Then, by considering the balance of the virtual works integrated in the space-time domain, the variational formulation in each space-time slab S n expressing the dynamic equilibrium and the displacementvelocity compatibility reads as: where V (S n ) is the space of kinematically admissible fields (w u , w v ), [·(t n−1 )] = (·)(t + n−1 ) − (·)(t − n−1 ) denotes the jump quantity in time at t n−1 and (·, ·) D denotes the inner product between two vector fields or two tensor fields integrated over a domain D, e.g.
It is worth recalling that to obtain (6), the integration by parts is also made in time, then a numerical flux between two successive space-time slabs should be chosen, as the fluxes are not continuous, e.g. at t n , we have (v(t − n ), σ(u(t − n ))) = (v(t + n ), σ(u(t + n ))). Due to causality, a natural choice is to take the numerical fluxes equal to the upwind fluxes (v(t − n ), σ(u(t − n ))). Now, let us consider the total energy functional E (u(t), v(t)) over the whole studied domain that is defined as follows: The following result concerning the stability of the time dG method (6) can be proved.
In the sense that (9) holds also for finite element solutions whatever the space-time mesh used to discretize each space-time slab S n , the time dG method (6) is unconditionally stable.
Hence, the kinetic and elastic energies dissipated in the displacement and velocity jumps in time between two successive space-time slabs guarantee the unconditional stability of the time dG method.
Otherwise, we recall that space-time meshes generally used for the time dG method are composed of a classical finite element mesh in space and one linear element in time in each space-time slab and so an implicit solver is actually obtained [4][5][6][7][8]17].

Velocity-stress space dG method
The governing equations of the same model problem of elastic wave propagation defined in the preceding are firstly given within the framework of the velocity-stress space dG method. The unified and wave oriented variational framework proposed in [17] is recalled. Then the exact upwind numerical fluxes developed in [17,18] are given before the presentation of the main result of the present work: the stability of the space dG method in the general case of anisotropic heterogeneous media with physical interfaces of discontinuous material properties.

Strong and variational formulations
We consider the same model problem of elastic wave propagation defined in the preceding section. However, within the framework of the velocity-stress space dG method proposed in [10,11,19] (see Eqs. (22.14-22.15) in [19]), its governing equations are put into the following strong first-order hyperbolic system with as primary unknowns the velocity and stress fields (v, σ): where the space derivation operator A ∂ x and its transpose (useful hereafter) are defined as follows: We note that in this section no body force is considered in the equilibrium equation without loss of generality of the purpose of the present work.
The tensorial compact form in (10) was proposed in [17]: the generalized unknown According to the definition of the space derivative operators (3), it is easy to show that, on the boundary ∂D of any subdomain D ⊆ , the flux operator F n (with n = n i e i the outward unit normal vector defined on ∂D) associated to the first-order system (10) is in fact equal to A n the Jacobian operator in the n direction, that is: In (13), the subscript index "n" indicates the dependency of F n and A n on n. In the following, the local orthonormal basis defined on ∂D will be denoted by (n, {t α } α=1,...,d−1 ). Furthermore, it is useful to recall one important result given in [17]: A n can be decomposed using its two eigenbases as follows: where longitudinal and quasi transverse wave modes propagating in the n direction. The subscript indices "qL" and "qT " respectively refer to terms "quasi longitudinal" and "quasi transverse". According to [17], the right and left eigenmodes corresponding to the nonzero eigenvalues of A n read as: with z ± n,k = ρλ ± n,k the acoustic impedance, w n,k = 1 ..,d−1 the unit eigenvectors of the usual eigensystem of the Christoffel tensor n : The definition of the Christoffel tensor n is recalled in the following: We recall that the word "quasi" means that in the general anisotropic case we have neither pure longitudinal wave mode verifying γ n,qL n nor pure transverse waves modes verifying γ n,qT ⊥ n, in contrast to the isotropic case.
To develop the variational formulation for the space dG method applied to the strong formulation (10), the main idea is to seek an approximated solution U h = (v h σ h ) T that are discontinuous in space, so we need to use a finite element mesh and the obtained space dG formulation is already discretized in space.
Let M h = { k } k denote a finite element mesh of . In the following, any element k of the mesh M h will be denoted by E and any of its neighboring elements by E . The discontinuous solutions in E and E are respectively denoted by U h and U h . Then the space dG variational formulation of the elastic wave model problem (10) for any elemnt E can be put into the following two equivalent forms: In (18),F n (U h , U h ) denotes a numerical flux used to replace the discontinuous flux F n (U h ) on the element boundary ∂E when the variational formulation (18a) is established. Hence, an appropriate choice of the numerical flux is essential for the success of the space dG method.
In the following, is firstly recalled the definition of numerical fluxes on the interior element boundary ∂E int = ∂E\(∂E ∩ ∂ ), which should also depend on the solution U h in the neighboring elements E of E. As for the exterior element boundary ∂E ext = ∂E∩∂ , the definition of the numerical fluxF n (U h ) should take into account the boundary conditions U h . It was given in [17] and will be recalled in Lemma 1.

Exact upwind numerical fluxes
Let us consider the interface of two adjacent elements E and E having respectively (ρ, C, U h ) and (ρ , C , U h ) as densities, elastic moduli and initial states (Fig. 1). The two outward unit normal vectors of E and E on their interface are respectively n and n and verify n + n = 0.
To define exact upwind numerical fluxes, we should consider the Riemann problem governing the states that are results of the propagation of the discontinuity U h − U h [15,[17][18][19][20]. In the following all the equations are written in the 3D case without loss of generality.
As discussed in [17], solving exactly the Riemann problem at element interfaces is not a sufficient condition to get physically sound numerical solutions. Indeed, different equivalent strong forms of the elastic wave problem give rise to different forms of the Riemann problem and consequently to different interface conditions. For example, the Riemann problem directly derived from the strong form (10) leads to the following interface conditions: Obviously, on a physical interface with ρ = ρ or/and n = n , (19) are generally not coherent with the classical interface conditions, i.e. the following velocity and stress vector continuities: In the present work, only the exact upwind numerical fluxes developed in [18] that are coherent with (20) are considered.
Recalling the results given in [18], the Rankine-Hugoniot jump conditions leading to the interface conditions (20) read as: It is worth noticing that the Riemann problem (21) corresponds to another equivalent form of the first-order velocity-stress system (10): Furthermore, by recalling (15) the definition of (R − n,k , L − n,k ) and (14) the decomposition of A n , it can be shown that: Solving of the Riemann problem (21) leads to the determination of the six unknown Then, the exact numerical fluxes are defined in the following way: We remark that (25b) gives another flux definition that is equivalent to (25a). It will be used in the proof of Lemma 1.
In [18], it has been proved that {α k ,α k } k=qL,qT 1 ,qT 2 are the solution of the following linear system of equations: In (26) respectively denote the harmonic and arithmetic means of the acoustic impedances of the k-th eigenvector. {˜L − n,k ,˜L − n ,k } are the perturbed left eigenmodes of {A n , A n } calculated by using the material properties of the adjacent element E in the following way: Two simpler but interesting cases to consider are: • Continuous case On an element interface with continuous material properties, it is straightforward that:α • Discontinuous isotropic case On an element interface separating two different isotropic materials, by simply recalling that in the isotropic case γ n,qL = n, γ n,qT 1 = t 1 and γ n,qT 2 = t 2 , it is straightforward that: We note that, only in the general case of an element interface separating two different anisotropic materials, we have

Stability analysis
Now the main result concerning the stability of the space dG method (18) using the exact upwind numerical fluxes (25)-(26) is given. Firstly, the stability of the space dG method is proved in both continuous case and discontinuous isotropic case whatever the space dimension. Then, in the general case of anisotropic media with physical interfaces separating different materials, the proof, which is much more complicated, is given only for the space dimension equal to 2.
Hereafter, EE = ∂E ∩ ∂E denotes the interface between any pair of adjacent elements E and E and the following notations are used for space jumps across it: Let us consider the total energy functional E (t) over the whole studied domain that is defined within the framework of the space dG method as follows: Firstly, the following lemma can be proved.

Lemma 1
By assuming that there is no source term inside ( i.e. f = 0) and the homogeneous Neumann and Dirichlet boundary conditions ( i.e. g = 0 and u D = 0), we have Proof Let us consider any element E. When ∂E ext = ∅, let us assume, as in [17] ("Exact upwind numerical fluxes" section), that the Dirichlet boundary conditions are strongly imposed by eliminating the prescribed degrees of freedom of velocity in the final system to solve. On the other hand, as also discussed in [17], it is natural to choose the following numerical flux on ∂ N where the Neumann boundary conditions are prescribed: Then, when the strategy of weak verification of the Neumann boundary conditions is chosen, the two equivalent space dG variational formulations (18) become: and, when the strategy of strong verification of the Neumann boundary conditions is chosen, i.e. σ h · n = g, they become: Now let us sum up the two equivalent space dG variational formulations (36) or (37), divide the obtained equation by 2 and take then, by remarking that the following equations are true: we finally get: Now, by summing up all the elements, remarking that: and taking into the definition (25b) of the numerical fluxesF n (U h , U h ) and (24a), the Eq. (34) can be finally proved.
Now, using Lemma 1, the following results concerning the stability of the space dG method can be proved.

Theorem 2
In the continuous case, the space dG method is stable in the sense that: with |Ã n | = |z ± n,k |L ± n,k ⊗ L ± n,k .
Proof In the continuous case, {α k ,α k } k=qL,qT 1 ,qT 2 are given by (30) and the following equations hold: Using the decomposition (24b) ofÃ n , the equation in (42) is straightforward. As |Ã n | is symmetric and positive-definite, the inequation in (42) is obvious.

Theorem 3
In the discontinuous isotropic case, the space dG method is stable in the sense that: Proof In this case, {α k ,α k } k=qL,qT 1 ,qT 2 are given in (31), then we get: with A * = |z − n,k |L − n,k ⊗˜L − n,k + |z − n ,k |L − n ,k ⊗˜L − n ,k . According to (15), (29) and (28), the definitions of (L − n,k , L − n ,k ), (˜L − n,k ,˜L − n ,k ) and (C − z,k , C − z,k ), we get: Herein, the orthonormal basis on EE chosen for E (resp. E ) is (n, {t l } l=1,d−1 ) (resp. (n , {t l } l=1,d−1 )) with t l + t l = 0. Then by taking into account the fact that in the isotropic case the eigenvector {γ n,k } k=1,d (resp. {γ n ,k } k=1,d ) can be taken equal to the elements of the corresponding basis, we have γ n ,k = −γ n,k (∀k). Putting the latter and (46) in (45), it is easy to get (44).
We note that a similar result for the velocity-strain space dG formulation was proved by Wilcox et al. in [15], similar roles played by the harmonic and arithmetic means of the acoustic impedances |z − n,k | R and |z − n,k | V were highlighted. It can be noticed that the proof of the stability in the discontinuous isotropic case is more tricky than the one in the continuous case. However, thanks to the fact that the relation γ n ,k = −γ n,k (∀k) remains true, the proof remains relatively simple, contrary to the stability analysis in the discontinuous anisotropic case. In the following, only the 2D case is considered and the stability is proved only under a sufficient condition.
In the discontinuous anisotropic case, Theorem 4 shows that the stability of the space dG method is much more complicated to establish, even in the 2D case. The difficulty is essentially due to the interactions between wave modes of different types, i.e. γ n,1 with γ n ,2 and γ n,2 with γ n ,1 , coming from two adjacent elements separated by a physical interface. The kinetic and elastic energies dissipated in the space jumps of velocity [v h ] and of stress vector [σ h · n] can still guarantee the stability, but under a sufficient condition. For the moment, it has not been possible either to remove this condition or prove that it is also a necessary condition. However, we remark that all equations developed in Theorem 4 and in its proof are perfectly symmetric with respect to all quantities coming from E and E , hence, we believe that it would be difficult to further improve the condition proposed in Theorem 4.
Concerning the sufficient condition (49), it depends obviously on the degree of discontinuity across a physical interface, as R 12 , V 12 and ζ are proportional to the latter. Future theoretical and numerical investigations are necessary to better understand this condition. However, a first analysis suggests that it should be satisfied for relatively high degrees of discontinuity.
In conclusion, for both time dG and space dG methods, the mechanism to ensure the stability is the energy dissipation of velocity, displacement or stress vector jumps in time or in space. However, contrary to the stability of the time dG method that is unconditional whatever the space and time discretizations used, the stability of the space dG method is proved for its variational formulation discretized in space but not yet discretized in time. When a time discretization scheme is applied, its own stability condition should also be taken into account. For example, in the present work, the standard four-stage fourth-order Runge-Kutta (RK) method is used for the time discretization of the space dG method, it is an explicit scheme submitted to Courant-Friedrichs-Lewy (CFL) type stability condition [17]. Inversely, even when the stability condition is satisfied by the time discretization method, if the stability is not proved at the continuous level for the space dG variational formulation, typically in the case of the use of approximated numerical fluxes, the stability of discretized solutions cannot be guaranteed.

Numerical investigations and comparisons
This section presents several numerical examples of comparison between both time and space dG methods. As previously indicated, concerning the time dG method, a space-time finite element mesh composed of a classical finite element mesh in space and one unique linear element in time is used to discretize each space-time slab. For the comparison, a same space finite element mesh is used by both methods, but the space dG method uses the standard four-stage fourth-order RK algorithm that is explicit and conditionally stable.
The strategy adopted in our works to choose space and time discretization parameters was detailed in [17]. In summary, the element size h E is determined by the shortest wavelength of interest, which is defined by the highest frequency of interest. Concerning the time step t, it is determined by the element size and the following CFL (Courant-Friedrichs-Levy) type conditions: for the space dG method it is a stability condition, while for the time dG method it is a condition to prevent high frequency modes of interest from numerical damping: In (52), N p is the order of FE basis function and CFL sDG the CFL number in classical sense. For the numerical examples presented hereafter, only linear (1D) or bilinear (2D, quadrilateral element with 4 nodes) elements are used, and we take CFL sDG ≈ C tDG−damping with N p = 1. So the time step used for the time dG method is three times larger than the one used for the space dG method.
As far as the external loadings are concerned, their dependency upon the time t is chosen to be a ricker signal, defined as follows: with T r the period of the ricker signal. We recall that the ricker signal (53) provides a perfectly controlled frequency content, its center frequency is f max = 2T −1 r and its cutoff frequency can be reasonably considered as f c = 3f max (see Figure 2 in [17]).

1D case with a physical interface
A 1D elastic rod = 1 ∪ 2 composed of two parts of different materials is studied. It was already considered in [17]. The density ρ = 2500 kg m −3 is uniform in the entire rod, while the Young's modulus in the two subdomains is respectively (from left to right) E 1 = 22.5 GPa and E 2 = a 2 E 1 . Hence, the phase velocity in the two subdomains is respectively equal to c L1 = 3000 m s −1 and c L2 = a c L1 . The total length of the rod is L = L 1 + L 2 = (1 + a)L 1 , with L 1 = 900 m. The impedance contrast between the two subdomains is It is obvious that the smaller the value of a, the higher the degree of discontinuity. Three degrees of discontinuity are considered, corresponding to three values of a: 1 (homogeneous case), 0.5 (discontinuous case) and 0.01 (highly discontinuous case).
A pressure loading g(t) is applied on the left end of the rod and the free boundary condition is prescribed on the right end. g(t) is a ricker signal with period T r = 0.04 s and cutoff frequency f c = 125 Hz. The corresponding shortest wavelength to consider is L (f c ) = 24 m. The time that waves need to propagate from the left end to the physical interface is equal to 0.2 s, so is the time from the physical interface to the right end. The total time of analysis is taken equal to 1.6 s, which corresponds to the time required for 2 round trips.
Firstly, the following element sizes h E 1 = 0.6 m and h E 2 = a h E 1 are used in the two subdomains, which results in N e = λ min h E = 40 elements in the shortest longitudinal wavelength of interest. The time step t = 62.5 µs is used by the space dG method and t = 20 µs by the time dG method, leading to CFL sDG ≈ C tDG−damping ≈ 1. Decreasing over time of the total energy in the entire rod is studied using two quantities: the first one is the total loss accumulated up to a time step: E(t)−max (E) max(E) (Fig. 2b), and the second one is the incremental loss between two successive time steps: Fig. 2b). Figure 2a leads to a first conclusion: with only one physical interface, there is no significant influence of the degree of discontinuity, i.e. of a, on the numerical damping rate over time. Moreover, with our choice of discretization parameters (h E , t), i.e. CFL sDG ≈ C tDG−damping ≈ 1, same levels of damping are observed for the two methods.
However, detailed analyses of the energy decreasing between two successive time steps (Figs. 2b, 3a) show that, at the physical interface, perturbations are observed only for the time dG method, which increase with decreasing value of a, while the space dG method shows a perfect behavior at the physical interface. Otherwise, due to the weak verification of the free boundary conditions, the energy damping of both methods is perturbed and such perturbations do not depend on a and are much more important for the space dG method than for the time dG method (Figs. 2b, 3b). It will be shown in the next section that such perturbations may become problematic in the 2D case. In conclusion, compared to the time dG method, the behavior of the space dG method across the physical interface is very interesting especially in cases involving a large number of physical interfaces. But future studies are necessary to understand and improve its behavior at external boundaries.
Secondly, influences of space-time discretization parameters (h E , t) are investigated. For both methods, the energy damping increases rapidly with mesh coarsening in space and in time (Fig. 4a). As expected, the energy damping of the space dG method depends principally on the element size in space (Fig. 4b), while the one of the time dG method depends principally on the time step (Fig. 4c).

2D case with physical interfaces
For the 2D case, the studied 2D rectangular elastic domain is presented in Fig. 5a and the width of the whole domain is L x = 5 mm. A similar but four times larger (in each space dimension) domain was already considered in [17]. The following four cases are simulated:  This case corresponds to a single-phase polycrystalline material with 85 elliptic grains of size 480 µm×240 µm (Fig. 5b). Each grain has the same Hooke tensor C given in the case (3)  Elastic waves are generated by applying a pressure loading g(x, t)e y in the vertical direction e y on an extremely short segment L e of length 50 µm at the center of the external top boundary of , while free surface boundary conditions are prescribed on the remaining parts (Fig. 5a). The pressure g(x, t)e y is uniformly distributed in space and its period of ricker signal in time is T r = 0.64 µs, resulting in f max = 3.1 MHz and f c = 7.8 MHz.
Bilinear square finite elements of size h E = 25 µm are used to discretize the whole domain . For the two first cases, this choice corresponds to respectively 12 and 18 elements in the shortest transverse wavelengths of interest T 1 (f c ) = 314 µm in 1 and  [21,22]. In particular, it has been shown that the mesh convergence depends not only on the wavelength to element size ratio but also on the grain size to element size ratio. The time step used for the space dG method is t sDG = 0.533 ns and the one for the time dG method is t tDG = 3 t sDG = 1.6 ns, leading to the value of CFL sDG = C tDG−damping equal to 0.37 for the case (1), to 0.56 for the case (2) and to 0.42 for the last two cases.
Evolution over time of the total, the kinetic and the elastic energies calculated by both space and time dG methods are compared in Figs. 6, 7. Both methods give very close curves for the kinetic energy, which is coherent with de comparison presented in [17] done only with the data of the velocity unknown field v h .
However, as the space dG method implemented in the present work uses the strategy of weak verification of the Neumann boundary conditions, the elastic energy calculated by the space dG method is highly perturbed by the non verification of those boundary conditions. Indeed, according to (40), for example, when the free Neumann boundary conditions are not exactly verified, the quantity (v h , g) ∂ N , with g = σ h · n, should be added in the energy evolution. In Fig. 6b, it is shown that this perturbation increases in the domain 2 with larger elastic moduli. In Fig. 7b, it is shown that this perturbation highly increases in the polycrystalline case, probably because the grains constantly reflect waves in the boundary layer and thus reinforces the perturbation. Future studies are therefore necessary to improve the stress field solution of the space dG method with respect to the Neumann boundary condition, by adopting for example the strategy of strong verification of them. Afterward, quantitative comparison of the damping behavior between both space and dG methods would be performed.

Conclusions
Comparisons and analyses of the time discontinuous space-time Galerkin method and the space discontinuous Galerkin method applied to the elastic wave propagation in anisotropic and heterogeneous media were presented. More particularly, the stability of both methods was investigated. It has been shown that the mechanism to ensure the stability of both methods is similar, in the sense that it is based on the energies dissipated in time or space discontinuities of unknown fields, but essentially different. Indeed, contrary to the stability of the time dG method that is unconditional whatever the space and time discretizations used, the stability of the space dG method is proved for its variational formulation discretized in space but not yet discretized in time. When a time discretization scheme is applied, its own stability condition should also be taken into account.
As the main contribution, a new result of stability of the velocity-stress space dG method was given. The stability of the space dG method was proved at the continuous level of time in the case of elastic media with continuous properties and the case of isotropic elastic media with discontinuous properties, whatever the space dimension. However, the proof of the stability was given only in the 2D case for anisotropic elastic media with discontinuous properties under a sufficient condition. Future theoretical and numerical investigations are necessary to better understand this condition. However, a first analysis suggests that it should be satisfied for relatively high degrees of discontinuity. Otherwise, the existence of such a sufficient condition in the 2D case suggests the existence of a similar condition in the 3D case. It should be important to develop it even if the proof in the 3D case would be very complicated.
Otherwise, the stability result of the space dG method was proved with the use of the exact numerical fluxes. It allows explaining previously observed phenomena [17]: when the degree of discontinuity at physical interfaces is high, the use of approximate numerical fluxes leads to an instability problem, which can be considered as intrinsic to the space dG method, in the sense that it cannot be resolved by simply reducing the time step used by the time discretization scheme.
Numerical applications to different elastic media with physical interfaces were presented and the numerical damping behavior of both space and time dG methods was investigated. In the 1D case, the quantitative comparison between the two methods was complete and has shown, at physical interfaces, the perfect behavior of the space dG method using the exact numerical fluxes. In the 2D case, the necessity to well satisfy the Neumann boundary conditions by the stress field calculated by the space dG method was highlighted. Future studies should be done to complete the quantitative comparison between the two dG methods and to assess the convergence rate of the space dG solver in terms of the energy and L2 norms and also of the physical quantities of interest, such as the attenuation and noise levels in the case of polycrystalline materials.
Then the proof of Theorem 4 is organized in eight parts.

• (2) Decomposition of ∂ t E into different parts to analyze
According to (26), we have: Using Lemma 1, we get: • (3) Calculation of the term concerning P EE According to (15) and (29) the definitions of (L − n,k , L − n ,k ) and (˜L − n,k ,˜L − n ,k ) and recalling that w n,k = 1 √ 2 γ n,k , we get:

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.