Je souhaiterais d'abord remercier les membres du jury d'avoir accepté de faire partie de mon jury. Je remercie Roland Martin et Stéphane Lanteri d'avoir accepté de rapporter sur ma thèse et d'y avoir apporté leur expertise scientifique et leurs remarques très pertinentes. Je remercie de plus Stéphane d'avoir accepté de présider le jury de cette thèse. Je remercie aussi Henri Calandra pour la dimension industrielle qu'il a apporté à ce sujet.

Je tiens également à remercier mes directeurs de thèse, Hélène Barucq et Christian Gout. Christian qui m'a suivi tout au long de mes trois années en tant qu'élève ingénieur à l'INSA de Rouen et puis durant mes trois années de thèse et grâce à qui j'ai pu avoir le stage ingénieur au sein de l'équipe Magique3D par lequel tout à commencer. Je remercie du fond du coeur Hélène qui a toujours été là quand j'en avais besoin et qui m'a permis de finaliser ces travaux. Pour tout ce qu'elle m'a apporté sur le plan humain je lui dis un grand MERCI. A ces deux personnes là je rajouterai Julien Diaz qui a dû supporter mes questions incessantes et qui m'a énormément apporté et m'a transformé en un développeur bien plus aguerri qu'à mon arrivé. Je remercie aussi l'ensemble de l'équipe Magique3D pour leur soutien indéfectible. Merci aux doctorants que j'ai côtoyé durant ces trois années. Les anciens Florian et ses demandes incessantes de tarte aux pommes, Elvira, Izar, Ajaral, Mamadou et les actuels Pierre, Rose, Vinduja, Nathan, Stéfano et la dernière, et non des moindres, Justine qui m'a supporté ces trois dernières années en tant que voisine de bureau et sans laquelle ces années n'aurait pas été les mêmes.

Le dernier paragraphe est dédié à ma famille et mes amis. Mes parents, qui m'ont toujours soutenu durant toutes ces années, mon frère et sa copine, mes grands-parents et tout ceux qui m'ont permit de devenir celui que je suis maintenant. Une pensée particulière à ma meilleure amie Marine qui a toujours été là pour me remonter le moral lorsque cela s'avérait nécessaire et qui m'a permis de surmonter toutes les épreuves que ces trois années ont engendré.

List of Figures

Introduction

Numerical seismic imaging is a topic that has evolved a lot thanks to the tremendous progresses of scientific computing together with the access to supercomputers. Seismic imaging begins with data acquisition which consists in recording reflected seismic waves generated by artificial sources and propagating in the underground. Then, receivers, which are usually placed at the surface, measure the direct waves and the reflected ones (see Fig- ure 1.1 for an illustration of an offshore survey campaign). To transform those data into Figure 1.1: Typical data acquisition for seismic marine imaging images of the underground, numerical techniques have been developed and distinguish themselves by the information they are able to retrieve and also by the degree of precision with which they provide with a reconstruction of the subsurface. From a numerical point of view, they are all engaged in a computationally intensive process that must be able to solve wave equations in large heterogeneous domains. Nowadays, two numerical methods have emerged for such a purpose. They are the Reverse Time Migration (RTM) and the Full Waveform Inversion (FWI). Both approaches are based upon the numerical solution of the full wave equation, which has to be performed many times (for each source and possibly the receivers), and the resolution is part of an iterative procedure that is supposed to converge towards the image of the domain under study. RTM is a method of migration by extrapolation of the full wave equation in the time domain while FWI is an inversion technique involving optimization methods adapted for solving ill-posed nonlinear inverse problems. RTM is a qualitative method of inversion that provides information on the location of reflectors. It is very popular because of its robustness and also because it is quite easy to implement but it only delivers information on the kinematic of the propagation medium. Hence it can be left for FWI which provides quantitative information by retrieving the dynamics of the domain through the reconstruction of its physical parameters.

This work aims at developing a propagator for acoustic and elastic waves that combines two different finite element methods. This new development is motivated by the objective of having a FWI in the time domain. FWI is implemented by following an iterative process which aims at updating the physical parameters of interest at each iteration (for instance the wave velocities) until convergence towards the image of the domain under exploration. A minimization method is carried out thanks to a gradient method applied to a misfit function that measures at each iteration the difference between measured and numerical data. The FWI has been introduced by Lailly [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations[END_REF]) and Tarantola (Tarantola, 1984a[START_REF] Tarantola | Linearized inversion of seismic reflection data[END_REF] in the time domain and it has been extended to the frequency domain by Pratt et al. [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF]. FWI is a numerical inversion method that is computationally intensive with still excessive assignments of the computational resources in realistic geophysical cases. Given that the method is sensitive to the accuracy of the numerical solution, it is problematic to give up accuracy in order to carry out the computations with available computing resources. Hence, it is of great interest of implementing a numerical method that involves less computational power than existing DGms. Magique-3D, in collaboration with Total, has developed a FWI software package in the 2D frequency domain [START_REF] Faucher | Contributions to seismic full waveform inversion for time harmonic wave equations: stability estimates, convergence analysis, numerical experiments involving large scale optimization algorithms[END_REF] that turns out to be less expensive than in the time domain because the images are carried out thanks to the inversion of a discrete number of frequencies. In 3D, one cannot conclude in that way because in the frequency domain, the numerical method is based on the use of linear solvers which are not yet able to address large problems like those considered for geophysical applications. There is thus an interest of developing FWI based upon time-dependent wave equations and this work aims at preparing the development of FWI in time domain by providing an accurate numerical method for simulating wave problems with affordable computational costs. In particular, because FWI implies solving several wave equations, the occupation of the memory must be controlled and for this reason, explicit time integration is preferred.

There exist several numerical methods that have demonstrated their efficiency for solving wave problems. In the geophysical community and in particular for industrial developments, the most common numerical approach is the Finite Difference method (FDm). It works with cartesian grids and it is well-known for providing numerical results with reasonable computational costs besides being quite easy to implement. However, regular grids are not able to easily reproduce the possible variations of the topography. For instance, local refinements can be carried out to fit with the topography variations but numerical dispersion can appear and pollute significantly the numerical solution. In such situation, the flexibility of Finite Element methods (FEm) with tetrahedral meshes, is a great advantage even if FEms require higher computational resources. This is true for finite element approximations based upon Lagrange polynomials. They have a mass matrix which is not diagonal with a stencil increasing with the order of approximation. In the case of an explicit time integration, FEm involves then the inversion of a matrix at each time step which increases significantly the computational costs and thus accentuates the interest for finite differences. Mass lumping could have been used but it turns out that mass lumping process is not able to keep the order of the Lagrange finite element method, as soon as its order exceeds two. To overcome this difficulty, Spectral Element methods have been investigated and Dimitri Komatitsch and his co-authors (see (Komatitsch & Vilotte, 1998;[START_REF] Komatitsch | The spectral element method for elastic wave equations: application to 2D and 3D seismic problems[END_REF]Komatitsch, Barnes, & Tromp, 2000a) among others) have clearly demonstrated that this method is very efficient for solving wave problems applied to seismology. This method distinguishes itself from other finite element methods by employing quadrature formulas which use Gauss-Lobatto-Legendre points. If the degrees of freedom of the method coincide with these quadrature points, the mass matrix is diagonal. This is an important feature which allows using explicit schemes for time integration. Then, it is worth noting that high order SEms must be developed on hexahedral grids, which tends to limit the use of this approach in an industrial context, at least until now, because the generation of meshes with hexahedra is difficult and the commercial facilities are rare. We can cite Hexamesh (https://www.dps-fr.com/maillage-hexa-automatique-catia), which has been developed recently, but it is devoted to applications that dispose of computer aided design and this is not the case of geophysical exploration that aims at retrieving the propagation domain itself. This has been a motivation for developing wave propagators which are based upon Discontinuous Galerkin (DG) approximations. DG methods (DGm) are now widely used for solving wave problems (see [START_REF] Hu | An analysis of the discontinuous Galerkin method for wave propagation problems[END_REF]Cockburn, Karniadakis, & Shu, 2000;[START_REF] Dumbser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes. ii. the three-dimensional isotropic case[END_REF][START_REF] Grote | Discontinuous Galerkin finite element method for the wave equation[END_REF][START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation[END_REF][START_REF] Bernacki | Time-domain parallel simulation of heterogeneous wave propagation on unstructured grids using explicit, nondiffusive, discontinuous Galerkin methods[END_REF]Käser, Dumbser, De La Puente, & Igel, 2007; de la Puente, Käser, Dumbser, & Igel, 2007; de la Puente, [START_REF] De La Puente | Discontinuous Galerkin methods for wave propagation in poroelastic media[END_REF][START_REF] Käser | Quantitative accuracy analysis of the discontinuous Galerkin method for seismic wave propagation[END_REF][START_REF] Chung | Optimal discontinuous galerkin methods for the acoustic wave equation in higher dimensions[END_REF][START_REF] Petersen | A space-time discontinuous galerkin method for the solution of the wave equation in the time domain[END_REF][START_REF] Delcourte | A high-order discontinuous Galerkin method for the seismic wave propagation[END_REF][START_REF] Wilcox | A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media[END_REF][START_REF] Dupuy | Wave propagation in heterogeneous porous media formulated in the frequency-space domain using a dis-References continuous Galerkin method[END_REF][START_REF] Hermann | Non-conforming hybrid meshes for efficient 2-d wave propagation using the discontinuous galerkin method[END_REF][START_REF] Boillot | Contributions à la modélisation mathématique et à l'algorithmique parallèle pour l'optimisation d'un propagateur d'ondes élastiques en milieu anisotrope[END_REF][START_REF] Faucher | Contributions to seismic full waveform inversion for time harmonic wave equations: stability estimates, convergence analysis, numerical experiments involving large scale optimization algorithms[END_REF] for instance). They are built on arbitrarily shaped meshes including tetrahedra and hexahedra. They consist in expanding the approximate solution in a basis of discontinuous functions that are defined locally in each element. By this way, whatever the order of the method is, one element communicates with its neighbors only and the mass matrix is block diagonal. This is an important feature of DGm since the mass matrix is easy to invert and the use of an explicit time scheme is possible with affordable computational costs. Moreover, the whole formulation takes advantage of the local definition of the basis functions by being implemented in parallel very naturally. As well, DG basis functions are locally defined and this favours the use of local orders of approximation. DGms are thus well suited for p and h adaptivity. In the framework of the partnership DIP between Inria and TOTAL, Magique-3d has developed propagators which use DGms, essentially for solving first-order wave equations. These methods have been tested in various situations and it turns out that they deliver very accurate solutions but their computational burden is clearly higher than the one of SEm (Komatitsch & Vilotte, 1998), even if most of the comparisons do not consider the same mesh for both methods. The results could thus be biased in the sense that even the grids do not contain the same number of nodes. Anyway, DGms provide very flexible approximation methods capable of handling propagation domains that include strong variations of the topography and are characterized by physical parameters that may vary inside the cells of the mesh. This is basically a situation in which DGm could perform better than SEm.

If we go back to Figure 1.1, we see that a typical example of geophysical exploration involves a propagation medium which is a layer of water over the subsurface composed of geological layers. It is clear that there is no interest in meshing the layer of water with tetrahedra. We can surely use a structured grid to compute the wavefield in this region.

The situation is different in the subsurface, where the geological layers may be arbitrarily shaped. Their discretization is obviously easier with tetrahedra and the idea of using an hybrid mesh becomes very interesting. However, employing a hybrid mesh does not address the question of developing a numerical method with reduced computational costs.

It is also necessary to elaborate a numerical method suitable for hybrid meshes. Given our knowledge of DGm and SEm, we propose here to couple these two numerical schemes.

It is worth noting that we have focused on the construction of a coupled method, which explains that we have not addressed the question of using general hybrid meshes and we acknowledge that this should be done in the future. Here we concentrate our efforts on the development and validation of a DG-SEm scheme for acoustic and elastic wave problems in time domain. It is worth noting that the interest for hybrid numerical methods is still increasing significantly because they offer a promising approach for tackling industrial problems. In that spirit, we can cite the work of N. Deymier [START_REF] Deymier | Etude d'une méthode d'éléments finis d'ordre élevé et son hybridation avec d'autres méthodes numériques pour la simulation électromagnétique instationnaire dans un contexte industriel[END_REF] who developed and implemented a hybrid method involving a DGm with a FDm (Yee scheme) for solving Maxwell equations in the time domain.

This manuscript is organized as follows. The next chapter provides a general setting of the wave equations solved for seismic applications. Given that the resulting propagator will be used for FWI, we will consider first-order wave equations because such formulations give access to the physical quantities of interest when dealing with inversion. In this chapter, we also describe the different types of media that can be encountered in geophysical exploration. Chapter 3 deals with the numerical schemes that are the basis of the numerical approach proposed in this work. First, we introduce two different time schemes that are the Leap-Frog and Runge-Kutta schemes. Then, we present the DGm in the case of elastodynamics and the SEm for acoustics. For each method, we introduce the variational formulation and the discretization space that gives rise to a semi-discrete system. We validate the implementation of each method in quadrangle meshes and this gives us the opportunity to compare DGm with SEm when both are applied in quadrangle meshes. We show that, in this case, SEm requires less computational resources but in the same time, DGm is more accurate. In Chapter 4, we develop a DG-SEm coupling scheme that is based on DG mixed-primal formulation allowing to handle the area restricted to the application of SEm as a DG macro-element. Chapter 4 is thus the heart of this work. For the sake of conciseness, we deal with the elastic wave equation in details and provide a short description of the acoustic case. In particular, by performing an analysis of the energy of the discrete system, we guarantee the stability of the numerical method when it is implemented with the Leap-Frog scheme. The method has also been implemented with the Runge-Kutta scheme of order 4 but, even if we did not observe any problem of stability, the corresponding energy study is still an ongoing work. We finish this chapter by addressing the case of the first-order elasto-acoustic equation that is employed for simulating the interaction between a fluid and a solid. For that purpose, the SEm is applied in the fluid area while the DGm is used in the solid area. It is worth noting that the interaction between the two different physics is modelled with a transmission condition that requires changing the formulation for achieving the coupling of SEm with DGm. In Chapter 5, we perform numerical experiments to validate the DG-SEm coupling scheme. For that purpose, we carry out comparisons with full DGm implemented with unstructured triangles in areas where it is more difficult to mesh the geometric and physical characteristics of the domain. All the experiments we have carried out clearly show that the DG-SEm computations are faster that the full DG ones with the same degree of accuracy. In the last chapter, we investigate the behavior of the numerical method when absorbing conditions are used. We address the case of absorbing boundary conditions (ABC) and the case of Perfectly Matched Layers (PML). Such boundary conditions have an essential role in the simulation of seismic waves since they allow performing regional computations which are less expensive than global simulations. Regarding ABCs, we introduce a first-order condition which has already been implemented in [START_REF] Boillot | Contributions à la modélisation mathématique et à l'algorithmique parallèle pour l'optimisation d'un propagateur d'ondes élastiques en milieu anisotrope[END_REF]. We extend this condition to the coupled formulation and illustrate its accuracy with some numerical experiments. More importantly, we address the question of stability regarding the use of PML. Actually, even in isotropic media where PMLs are known to be stable, we have observed that their integration in the full DG scheme can lead to instabilities. On the contrary, there are no surprises when using SEm. This gives us the idea of coupling DGm with SEm employed in the layer in which the PML acts. We then observe that such a coupling is successful in stabilizing the PML when using a DG scheme in the rest of the computational domain and for that purpose, we have implemented the ADE-PML [START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML)[END_REF]. This PhD thesis is based upon an important contribution to software development. Indeed, before its beginning, some fundamental bricks were missing in the pieces of software currently maintained in the team Magique-3D. For instance, the SEm has been developed from scratch and DGm was only implemented in triangular/tetrahedral meshes. The work has been done with twofold objectives: (1) implementation of a DG-SE method for solving wave equations in heterogeneous media; (2) validation of the method in view of disposing of a proof of concept in the context of seismic imaging. The development of a DG-SE method has been launched as an exploratory topic and we have decided to work first in the framework of Elasticus which is a software package devoted to the implementation of new methods requiring new data structures and this is particularly the case of DG-SE method. Hence as a preamble of the manuscript, we describe the main achievements that have been carried out regarding software development.

Software development

The numerical simulations displayed in this work have been performed with the software Elasticus, which is developed in the team Magique3D. At the beginning of this work, Elasticus was able to solve numerically wave equations using DGm in unstructured triangle meshes in 2D and unstructured tetrahedral meshes in 3D. Elasticus proposed three different propagation media: fluid domains in which acoustic waves propagate, solid domains in which elastic waves travel and the elasto-acoustic coupling. Regarding time discretization, Elasticus offers three different schemes: second-order Leap-Frog schemes (only for fluid and solid domains), Runge-Kutta of order 2 and 4 (for the three types of media). As far as boundary conditions are concerned, Elasticus is endowed with first order Absorbing Boundary Conditions (ABC) and Perfectly Matched Layers (PML) as well, using Auxiliary Differential Equation (ADE)-PML with Convolutional-PML or PML formulation [START_REF] Gedney | An auxiliary differential equation formulation for the complex-frequency shifted PML[END_REF][START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML)[END_REF]. These boundary conditions can be applied to the elastic wave equation and to the acoustic wave equation as well. In the framework of this thesis, several important features have been added to Elasticus. Since the beginning, we knew that we had to include SEm in Elasticus and that a comparison with DGm was necessary. Given that SEm works with quadrangles and hexahedra, we decided first to implement a DGm working with meshes composed of quadrangular cells in 2D and cuboids in 3D and with hybrid meshes composed of both triangles/tetrahedra and quadrangles/cuboids. Once these new functions have been added to Elasticus, we have implemented the SEm in 2D and 3D both for acoustic and elastic domains. It is worth noting that, at first, the development of SEm in elastic domains was not a priority for us since our objective was to apply SEm in fluid domains. Nevertheless, having a SEm for elastic domains was very interesting to go further in comparisons with DGm. Moreover, we have realized that it could be used to stabilize the PML in Elasticus, as we explain in Chapter 6.

Having full SEm and DGm, the main objective of this thesis has been addressed: we have constructed a numerical scheme including the coupling of SE and DG elements working with meshes composed of different areas having structured and unstructured cells. It is worth noting that structured cells are quadrangles in 2D and cuboids in 3D while unstructured cells are triangles in 2D, tetrahedra in 3D. This feature is justified by the fact that the coupling has been carried out in order to use DGm in the unstructured triangular/tetrahedral part and SEm in the quadrangular/cuboidal one. The coupling has been achieved to work with all the different time schemes that are available in Elasticus. Even if our main objective was the solution of elasto-acoustic problems with SE-DGm, the coupled scheme works in full acoustic and full elastic domains and for these two configurations, we can apply a Leap-Frog scheme and Runge-Kutta schemes while for solving the elasto-acoustic problem, we only use Runge-Kutta integrations required because the space discretization involves upwind fluxes. In terms of boundary conditions, we have implemented an Absorbing Boundary Condition both for DGm and SEm. It is worth noting that the ABC is plugged into the Sem and the DG-SEm formulations weakly, handling in the formulation as a penalization. In acoustics, we have used the first order ABC designed by Engquist and Majda in [START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF][START_REF] Engquist | Radiation boundary conditions for acoustic and elastic wave calculations[END_REF]. For elastic wave equation, we have used the first order ABC designed by [START_REF] Cohen | Computational methods for transient analysis[END_REF] for isotropic media, and the one derived in [START_REF] Boillot | Contributions à la modélisation mathématique et à l'algorithmique parallèle pour l'optimisation d'un propagateur d'ondes élastiques en milieu anisotrope[END_REF][START_REF] Barucq | Absorbing boundary conditions for 2D tilted transverse isotropic elastic media[END_REF] for transversely isotropic media. For the DGm/SEm coupling, we have taken advantage of the fact that we could consider the area in which the SEm is applied as a DG macro-cell (as we explain in Chapter 4). This approach eases significantly the coupling of the two ABCs. Finally, we have extended the Perfectly Matched Layer formulation implemented first for DGm and SEm apart to DG-SEm (see chapter 6 for details). All the above developments have been achieved both in 2D and in 3D.

Chapter 2

General setting of wave equations

In this chapter, we introduce the equations governing the wave propagation in a fluid (acoustic waves) and a solid (elastic waves) and how they are coupled in the case of propagation domains allowing the coupling between acoustic and elastic waves as depicted for instance in the introduction. We also describe the implementation of the source term. It is worth noting that we do not address in details the boundary conditions that we use for simulation. This is due to the content of Chapter 6 which focuses on boundary conditions.

Acoustic system: equation of propagation

The acoustic wave equation results from three fundamental laws represented by the Euler equation, the principle of mass conservation and the thermodynamics of the fluid. These equations link the pressure P , the velocity v and the density ρ which are the parameters characterizing the fluid. The Euler equation is derived as an application of the fundamental principle of the dynamics. In absence of external sources, it reads:

ρ(x, t) (∂ t v(x, t) + (v(x, t) • ∇)v(x, t)) = -∇P (x, t), (2.1)
where ∇ is the gradient operator. The equation governing the mass conservation is given by:

∂ t ρ + ∇ • (ρv(x, t)) = 0, (2.2)
where ∇• is the divergence operator. We now assume that the fluid is initially immobile, that is ρ(x, 0) = ρ 0 , P (x, 0) = p 0 and v(x, 0) = 0. Then, the fields can be parametrized as follows:

ρ(x, t) = ρ 0 (x) + δρ(x, t), P (x, t) = p 0 (x) + p(x, t). (2.3)
It is worth noting that p 0 is the static pressure which corresponds more generally to the pressure in absence of acoustic wave. Assuming that the acoustic perturbations are small, we can linearize the velocity and the perturbation of the density up to the first order and we get:

ρ 0 ∂ t v(x, t) = -∇p(x, t) (2.4) and ∂ t (δρ) + ∇ • (ρ 0 v(x, t)) = 0.
(2.5)

These two equations are not sufficient to fully characterize the fluid. We miss an equation of state linking the pressure and the density of the fluid. In the case of acoustic perturbations, the fluid is adiabatic and we have:

δρ = ρ 0 χp, (2.6)
where χ is the adiabatic compressibility coefficient which is a constant under hypothesis of small perturbations. Hence we can remove δρ from (2.4) to end up with the equation:

∂ t p + 1 χ ∇ • (v(x, t)) = 0.
(2.7)

Finally we obtain the linear acoustic system:

ρ 0 ∂ t v(x, t) = -∇p(v, t), χ∂ t p(x, t) + ∇ • v(x, t) = 0.
Most of the time, seismic imaging aims at reconstructing velocity models. That is why we prefer to rewrite the previous system by introducing the propagation velocity c 0 = 1 √ χρ 0 of acoustic waves. Hence, we will consider the first order velocity-pressure formulation of the acoustic wave equation:

         ρ 0 ∂ t v(x, t) = -∇p(x, t), 1 c 2 0 ρ 0 ∂ t p(x, t) + ∇ • v(x, t) = 0.
(2.8)

To be solved, this system should be completed with initial conditions and boundary conditions depending on the case we are considering. This question will be clarify in the following chapters.

Elastic waves: equation of propagation

The elastodynamic system is obtained from equations elaborated from the theory of strengths of materials combined with the equation of motion. There are two different formulations depending on whether the displacement or the velocity is considered as unknown. Regarding geophysical applications we deal with, we focus on the formulation using the velocity v.

Consider the propagation domain is being deformed from a given force. Then, it is elastic if it returns to its original shape when the effect of the force stops. In the case of linear elasticity, the deformation is supposed to be small. Then the strain tensor ǫ reads:

ǫ(v) = 1 2 ∇v + (∇v) T .
(2.9)

There are constraints which work to return the solid in its original shape. If the deformations are small, the strain tensor is connected to the stress tensor σ by the constitutive law:

∂ t σ(x, t) = C(x)ǫ(v(x, t)), (2.10)
where C denotes the elasticity tensor. As far as the equations are written component by component, we have

∂ t σ ij = k,l
C ijkl ǫ kl . The motion equation is deduced from the Newton law in absence of source term: .11) where ρ denotes the density. Then, by injecting (2.9) and (2.10) in the time derivative of (2.11), we end up with the elastodynamic system:

ρ(x)∂ t v(x, t) = ∇ • σ(x, t). ( 2
ρ(x)∂ 2 t v(x, t) = ∇ • C(x)ǫ(v(x, t) .
(2.12)

Therefore, we get a second order partial differential system with v as unknown. By considering the stress tensor σ as unknown as well, the second order equation transforms into a first order system of the form:

ρ(x)∂ t v(x, t) = ∇ • σ(x, t), ∂ t σ(x, t) = C(x)ǫ(v(x, t)). (2.13)
This is actually the formulation we solve and for which we aim at developing a hybrid numerical software package suitable for seismic imaging.

Principal seismic waves

Into the Earth, elastic waves can be P and S waves with possible conversions of P into S and vice versa. It is possible to characterize them from wave equation solutions in isotropic medium and in absence of source term. Let us introduce the scalar φ and the vector field Ψ as solutions of the wave equations

ρ ∂ 2 φ ∂t 2 -(λ + 2µ)∆φ = 0 (2.14) and ρ ∂ 2 Ψ ∂t 2 -µ∆Ψ = 0, (2.15)
where ρ is the density and λ, µ are the Lamé parameters (see Section 2.4 of this chapter.)

The solution to (2.14) is a scalar wave propagating with the velocity V P = λ + 2µ ρ while the solution to (2.15) is a wave vector field propagating with the velocity V S = µ ρ .

Let us now consider the second-order equation:

ρ ∂ 2 u ∂t 2 -µ∆u -(λ + µ)∇(∇ • u) = 0,
where u denotes the displacement field and ∇ is the gradient operator. It is easy to check that a solution u can be written through the Helmholtz decomposition

u = ∇φ + curlΨ,
where φ is solution to (2.14) and Ψ is solution to (2.15). Such a property mainly results from the well-known properties:

• the vectorial Laplace operator can be decomposed into ∆ = ∇∇ • -curlcurl;

• curl∇ = 0;

• ∇ • curl = 0.
Then, the component u P = ∇φ corresponds to a longitudinal wave or Primary (P) wave as the one that propagates faster and is curl free. We can observe that the component u S = curlΨ is a transverse wave or Secondary (S) wave which is divergence free. S waves propagate more slowly than P waves.

Characterization of elastic media

In this work, we address both isotropic and anisotropic media. Geophysical media are generally anisotropic with possible isotropic areas. More precisely, we will consider three types of elastic media: isotropic, Vertical Transverse Isotropic (VTI) and Titled Transverse Isotropy (TTI). It is the elasticity tensor C that contains information about the properties of the medium under study. It is a tensor of order four whose components are denoted by C = (C ijkl ) 1≤i,j,k,l≤d where d is the dimension of the domain. According to Hooke's law, it is defined from the deformation tensor ǫ and the stress tensor σ. In the most general case, it has 81 coefficients. It is symmetric because ǫ and σ are symmetric. Moreover, assuming that σ derives from a potential, the elastic tensor becomes invariant under permutations of indices. This last property allows us to reduce its number of coefficients to 21. Hence, using the Voigt notation where the indices are gathered by pairs and replaced following a correspondence table, we can rewrite C as a matrix simply denoted by C which is symmetric and has at the most 21 independent entries. Positiveness is also an important property for the elasticity tensor, that is:

∀ symmetrical tensor ξ, (Cξ) : ξ ≥ 0
where ":" is the tensorial scalar product.

Isotropic medium

In the case of an isotropic medium, the elasticity tensor C contains only two independent coefficients known as the Lamé constants (λ, µ), which can be written in terms of Poisson coefficient ν and Young modulus E as follows:

λ = Eν (1 + ν)(1 -2ν) , µ = E 2(1 + ν) .
(2.16)

The Young modulus E measures the stiffness of a solid material and the coefficient of Poisson ν gives a measurement of the expansion of the material in the directions perpendicular to the direction of compression. The Lamé coefficients allow to express the stress tensor σ as a function of the deformation tensor ǫ:

σ = 2µǫ + λtr(ǫ)I,
where tr stands for the trace of matrices (that is the sum of the diagonal entries) and I is the identity tensor.

In an isotropic medium, the matrix C has then the following form:

C iso =          λ + 2µ λ λ 0 0 0 λ + 2µ λ 0 0 0 λ + 2µ 0 0 0 µ 0 0 µ 0 µ         
.

(2.17)

An isotropic medium can also be parametrized by considering the propagation velocities of P and S waves denoted respectively by V p and V s and the density of the medium ρ. In that case, the matrix C reads:

C iso =          C 11 C 11 -2C 66 C 11 -2C 66 0 0 0 C 11 C 11 -2C 66 0 0 0 C 11 0 0 0 C 66 0 0 C 66 0 C 66          , ( 2.18)
with C 11 = ρV 2 p and C 66 = ρV 2 s .

Vertical Transverse Isotropic medium

In realistic geological media, the waves do not generally propagate uniformly in every direction. For instance, the propagation domains include orthotropic materials whose properties differ along three mutually-orthogonal twofold axes of rotational symmetry. Such cases provide examples of anisotropic media whose properties change when measured in different directions. In the particular case where the axis of the frame of study and the plan of orthotropy are colinear, the matrix C is represented by 9 independent coefficients:

C ortho =          C 11 C 12 C 13 0 0 0 C 22 C 23 0 0 0 C 33 0 0 0 C 44 0 0 C 55 0 C 66          . (2.19)
In seismic imaging, we consider the polar anisotropy also called Transverse Isotropy (TI). It means that orthotropy only occurs in one direction with only one symmetrical axis usually very close to the vertical. In this case, the wave front is isotropic in a plan orthogonal to this axis. Figure 2.4 depicts the difference between an isotropic, a TI and an orthotropic medium. In an orthogonal frame where e z coincides with the axis of symmetry (denoted by ẽz ), the polar anisotropy is called Vertical Transverse Isotropy (VTI) and the elasticity matrix is written with 5 independent constants:

C V T I =          C 11 C 11 -2C 66 C 13 0 0 0 C 11 C 13 0 0 0 C 33 0 0 0 C 44 0 0 C 44 0 C 66          , ( 2.20)
where:

                   C 11 = ρV 2 p (1 + 2ǫ), C 33 = ρV 2 p , C 13 = ρ (V 2 p -V 2 s ) 2 + 2V 2 p δ(V 2 p -V 2 s ) -ρV 2 s , C 44 = ρV 2 s , C 66 = ρV 2 s (1 + 2γ).
(2.21)

The three constants ǫ, δ and γ have been introduced by Thomsen in [START_REF] Thomsen | Weak elastic anisotropy[END_REF]. They characterize a transverse isotropic medium and are equal to:

                           ǫ = C 11 -C 33 2C 33 , δ = (C 13 + C 44 ) 2 -(C 33 -C 44 ) 2 2C 33 (C 33 -C 44 ) , γ = C 66 -C 44 2C 44 .
(2.22)

If the three Thomsen constants are equal to zero, we retrieve C iso .

Tilted Transverse Isotropy (TTI)

In general, the orthogonal frame does not coincide with the frame of orthotropy. In that case, the transverse isotropy is tilted and the corresponding property of the medium is called Tilted Transverse Isotropy. The tilt is defined by two angles, θ and φ given by: θ = ( e z , ẽz ), φ = (e x , P xy ẽz ),

where P xy denotes the projection on the plane (e x , e y ) (see Figure 2.5). One interesting property of TTI is that it can be considered like a rotation of VTI. However, because the matrices of elasticity are derived from tensors, the rotation is represented in the tensorial algebra and it leads to two possible representations of the TTI tensor: either the tilt is written using the angles as a crossing matrix M (θ,φ) linking the VTI tensor and the TTI tensor, or the tilt is represented by using a rotation matrix R (θ,φ) .

In this work, we have adopted the second option with R (θ,φ) given by :

R (θ,φ) =    cos θ cos φ cos θ sin φ -sin θ -sin φ cos φ 0 sin θ cos φ sin θ sin φ cos θ    .
(2.23)

The TTI elasticity tensor is denoted by C and its coefficients are obtained from the ones of the VTI elasticity tensor C as follows:

Cijkl = d p=1 d q=1 d r=1 d s=1 R pi R qj R rk R sl Cpqrs 1 ≤ i, j, k, l ≤ d.
The resulting TTI matrix is composed of 21 independent coefficients written from 7 constants: V p , V s , ρ and the three Thomsen constants ǫ, δ and γ. However, even if they depend only on seven constants, the coefficients of the tensor do not have a simple 2.5. ELASTO-ACOUSTIC SYSTEM analytical formula, so by sake of simplicity, the TTI elasticity matrix is written as:

C T T I =          C 11 C 12 C 13 C 14 C 15 C 16 C 22 C 23 C 24 C 25 C 26 C 33 C 34 C 35 C 36 C 44 C 45 C 46 C 55 C 56 C 66         
.

Elasto-acoustic system

In the case of seismic marine (see Figure 2.6), the propagation domain is basically composed of a layer of water representing the ocean and a layer of solid materials standing for a piece of the sea bed. The typical propagation area is thus a square domain in two dimensions (or a cuboid in three dimensions) denoted by Ω. We note x ∈ R d (with d = 2, 3) the space variable and t ≥ 0 the time variable with t ∈ [0, T ]. Inside the sea, as formerly discussed in the introduction of the manuscript, we know that only P waves propagate while there are P and S waves which can travel inside the ocean floor. The P waves are solutions to the acoustic wave equation that we first introduce.

Figure 2.6: Typical data acquisition for seismic marine imaging

Implementation of source terms

In the framework of seismic exploration, we consider the propagation of waves which are generated by artificial sources. In the presence of a source, the source term S is taken into account in one of the two equations and in general, we consider source terms which vary in time and space independently. If the source term is included through the equation of motion, we have S(x, t) = s(x)s(t)

where s is a vector with the same dimension than v and s is a scalar. Then we solve the system

   ρ(x)∂ t v(x, t) = ∇ • σ(x, t) + s(x)s(t), ∂ t σ(x, t) = C(x)ǫ(v(x, t)) .
(2.24)

If the source term is included through the second equation of elastodynamics, we have

S(x, t) = s(x)s(t)
where s is a tensor and s is a vector with the same dimension than v. In that case, the system to be solved is:

   ρ(x)∂ t v(x, t) = ∇ • σ(x, t), ∂ t σ(x, t) = C(x)ǫ(v(x, t)) + s(x)s(t).
(2.25)

It is worth noting that in the case of an isotropic medium, the diagonal terms of tensor s generate P waves while the extra-terms create S waves.

For the numerical experiments, we use to have a time-dependent term which is a Ricker function. It is the second derivative of a Gaussian and depends on a peak frequency f peak . Its variation follows the formula:

s(t) = (1 -2π 2 f 2 peak (t -t peak ) 2 )e -π 2 f 2 peak (t-t peak ) 2 .
where t peak is the delay in time, usually set as 1.2/f peak . An example of Ricker source is depicted in Figure 2.7.

Boundary conditions

There are two classical boundary conditions which represent either a free surface or a rigid surface. The condition of free surface says that the traction vanishes, that is σn = 0, 2.7. BOUNDARY CONDITIONS n denoting the normal vector associated with the surface. It is the condition that is commonly employed for modeling the physical interface between air and a solid. The condition of rigid interface describes the interface between a solid and a rigid material. It is given by u = 0 on the surface where u denotes the displacement field. In the framework of this work, we restrict our use to the free surface condition. In acoustics, the free surface condition amounts imposing the pressure p is equal to zero on the boundary of interest.

In this work, we will consider media composed of fluid and solid areas. At the fluid-solid interface, it is then necessary to define adequate conditions of transmission. If Γ denotes the physical interface between the fluid and solid domain, the velocity v s in the solid domain is linked to the velocity v f and the traction is related to the pressure in the fluid domain as follows:

v f • n = v s • n (Γ), σn = -pn (Γ).
(2.26)

where n denotes the unitary normal vector defined on Γ and supposed to be oriented from the fluid area to the solid one. The natural conditions of transmission are thus the continuity of the normal speed and the continuity of the normal constraint.

To perform numerical simulations, it is necessary to truncate the computational domain.

The truncation boundary must be as transparent as possible to avoid the pollution of the numerical solution by reflected waves that it could have generated. Here, we use two different boundary conditions which are Absorbing Boundary Condition (ABC) and Perfectly Matched Layers (PML). The implementation of these boundary conditions along with their properties are detailed in chapter 6 of this manuscript.

Chapter 3

Key numerical methods

In this chapter, we introduce the numerical approximations that will be the key pieces of the numerical method developed in this work. We consider numerical integration both in space and time. The spatial discretization is based on finite elements and for the time integration, we favour explicit schemes of order 2 and 4.

Time discretization

The software package Elasticus offers the possibility of using three different time schemes which are Leap-Frog, Runge-Kutta of order 2 and Runge-Kutta of order 4. In this section, we provide a brief description of these schemes applied to a semi-discrete system obtained with a finite element approximation in space. Both the semi-discrete acoustic and the elastic wave equations can be written as:

d dt Y (t) = F (t, Y (t)) (3.1)
where Y is a pair of vectors (W, Q) corresponding to the finite element approximations of each unknown (velocity and pressure in acoustics, velocity and strain tensor in elastodynamics).

Leap-Frog scheme

The equation (3.1) can be written in the form

           d dt W = A 1 Q + A 2 W d dt Q = B 1 W + B 2 Q, (3.2)
where A 1 , A 2 , B 1 and B 2 are four matrices. The integration is carried out in the interval [0, T ] which is divided into time steps ∆t. Let W n be the approximation of W evaluated at discrete time t=n∆t, n ∈ N. We denote by Q n+ 1 2 the approximation of Q at the discrete time t = (n + 1 2 )∆t. Then, the Leap-Frog scheme consists in computing the next iterates in time by applying the following scheme:

             W n+1 -W n ∆t = A 1 Q n+ 1 2 + A 2 W n+1 + W n 2 Q n+ 3 2 -Q n+ 1 2 ∆t = B 1 W n+1 + B 2 Q n+ 3 2 + Q n+ 1 2 2 (3.3)
When A 2 and B 2 are zero, we obviously get an explicit representation of the discrete unknown since we can write W n+1 directly from the previous iterates W n and Q n+ 1 2 . The same accounts for Q n+ 3 2 which is given in terms of Q n+ 1 2 and W n+1 . This is the case when solving the acoustic wave equation with a spectral element method in space and Dirichlet or Neumann boundary condition on the external boundary. It is also the case with first order absorbing boundary conditions whose implementation leads to A 2 and B 2 diagonal. In the general case, the fact that the scheme is explicit depends on the algebraic structure of the matrices A 2 and B 2 which is determined by the underlying finite element discretization. For instance, when using a DG approximation involving centered fluxes, both A 2 and B 2 vanish and the scheme admits an explicit representation. But when upwind fluxes are required as for instance in the fluid-solid coupling, the scheme is implicit. As previously mentioned, the applications we are addressing require using explicit time integrations. This is due to memory limitations and justify that when applying a Leap-Frog scheme, we routinely employ a DG formulation with centered fluxes.

Runge-Kutta schemes

The Runge-Kutta (RK) scheme is based on the general form (3.1) of the semi-discrete system. It is based on the integration formula:

Y ((n + 1)∆t) -Y (n∆t) = (n+1)∆t n∆t F (t, Y (t))dt (3.4)
The order of RK scheme depends on the approximation order of the right-hand side of (3.4) which involves a quadrature formula. For that purpose, we apply first the following change of variable: t = (n + τ )∆t, where τ ∈ [0, 1]. Then we have:

(n+1)∆t n∆t F (t, Y (t))dt = ∆t 1 0 F ((n + τ )∆t, Y (n + τ )∆t))dτ (3.5)
Then, the integral along [0, 1] is approximated thanks to a quadrature formula. We use a quadrature formula with s points and we denote by τ i the quadrature points and ω i the quadrature weights. Then we have : (3.6) and:

Y ((n + 1)∆t) ≈ Y (n∆t) + ∆t s i=1 ω i F ((n + τ i )∆t, Y ((n + τ i )∆t)),
Y ((n + τ i )∆t) = Y (n∆t) + ∆t τ i 0 F ((n + τ )∆t, Y ((n + τ )∆t))dτ.
(3.7)

TIME DISCRETIZATION

In Equation (3.7), we interpolate F using still the points τ i as interpolation points. We then have:

τ i 0 F ((n + τ )∆t, Y ((n + τ )∆t))dτ = τ i 0 ψ j (τ ) s j=1 F ((n + τ j )∆t, Y ((n + τ j )∆t))dτ (3.8)
where ψ j , j = 1, ..., s are the Lagrange polynomials based on the points τ j , j = 1, ..., s. We finally get: In this work, we only consider explicit RK scheme and in particular Runge-Kutta of order two (RK2) and four (RK4) which have been implemented in Elasticus.

Y ((n + τ i )∆t) = Y (n∆t) + ∆t s j=1 a ij F ((n + τ j )∆t, Y ((n + τ j )∆t)) (3.9) with a ij = τ i 0 ψ j (τ )dτ. ( 3
RK2 scheme is constructed by using the midpoint quadrature formula which is given by

1 0 G(c)dc ≈ G( 1 2
) for any continuous function G ∈ [0, 1]. Then, we have:

Y ((n + 1)∆t) -Y (n∆t) = ∆t 1 0 F (t(n + τ )∆, Y ((n + τ )∆t))dτ ≈ ∆tK 2 , (3.11) with K 2 = F ((n + 1 2 )∆t, y((n + 1 2 )∆t)).
The RK2 is then given by:

K 1 = F (n∆t, Y (n∆t)) and
(3.12)

Y n+1 = Y n + ∆tK 2 with (3.13) K 2 = F ((n + 1 2 )∆t, y((n + 1 2 )∆tK 1 )). (3.14)
The Butcher table of RK2 is given by: 0 0.5 0.5 0 1

Table 3.3: Butcher table for RK2 scheme

In the same way, one can define the RK4 scheme which is represented by the following Butcher table: 0 0.5 0.5 0.5 0 0.5 1 0 0 1

1 6 1 3 1 3 1 6
Table 3.4: Butcher table for RK4 scheme

As for the Leap-Frog time-scheme, Runge-Kutta scheme stability is conditioned by a CFL number.

Discontinuous Galerkin formulation

In the team Magique-3D, we mainly use Discontinuous Galerkin (DG) approximations that have proved to be very efficient for reproducing the propagation of waves in heterogeneous media. We begin this section with an overview of DGm properties that led us to expend our efforts on this technology for solving seismic problems.

Overall setting of DGm

DG methods or DGms appeared in the 70's. They distinguish themselves from continuous Galerkin methods by the use of discontinuous basis functions, the discontinuity being at the interface between two elements. In DG approaches, each element has its own basis function set, which gives the possibility of working element by element. Actually, one element communicates with its neighbours only which makes the method well-suited for parallel computations. DGms are adapted to parallel computing in particular when an explicit time integration scheme is applied. Their efficiency is clearly significant when the interpolation order is high which makes DGms good for accurately approximating highorder partial derivative operators. Regarding the structure of the semi-discrete system, it is characterized by a block-diagonal mass matrix, hence with fast inversion, which eases the use of explicit schemes. It is particularly interesting for seismic applications that demand to perform a lot of simulations with a limited use of the memory normally reserved for inversion algorithm. Indeed, we have to keep in mind that we aim at elaborating a piece of software adapted to an implementation inside an inverse problem algorithm. Hence our main motivation is more reducing the computational cost than improving the accuracy of existing schemes. DGms can work with structured and unstructured meshes with various support mesh elements (hexahedra, tetrahedra, quadrangle, triangle) with conforming or nonconforming grids. This is an interesting feature regarding possible topography details that have to be considered when dealing with geophysical domains. The fact that the basis functions are discontinuous simplifies the implementation of hpadaptivity which opens up the possibility of using non-conforming meshes (h-adaptivity) and having different orders of approximation in each element (p-adaptivity). The padaptivity is quite simple to implement because it only requires to have a structure to store each order we use. It preserves the structure of the stiffness matrix and only modifies the size of the blocks inside the matrix. On the other hand, h-adaptivity implementation is more tricky. Indeed h-adaptivity requires changing the structure of the mesh by possibly introducing the so-called hanging nodes. A hanging node is a vertex of one element that also belongs to the interior of the edge of another element. Figure 3.2 shows an example of non conforming mesh with a hanging node in the middle of the edge of one element. This implies that it is impossible to use reference elements to simplify the calculation of the integrals defining the coefficients of the discrete system and by this way, the construction of the discrete system becomes more complex. The 3D case brings more difficulty because the intersection between two elements can take different forms, like a triangle or a quadrilateral. mesh with P 1 andP 2 At first, DGm has been used for solving elliptic and parabolic equations by Reed and Hill (see [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF]). DGm has been introduced for hyperbolic equations in the work of [START_REF] Lasaint | On a finite element method for solving the neutron transport equation[END_REF] and more developed ten years later in the ones of Nävert and Pitkäranta, [START_REF] Johnson | Finite element methods for linear hyperbolic problems[END_REF][START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF] where they established an optimal order of convergence for the method. A nice state of the art of DGm until the year 2000 can be found in [START_REF] Cockburn | The development of discontinuous galerkin methods[END_REF]. During the last twenty years, DGm became very popular for linear and non-linear hyperbolic equations. In particular we can refer to the work of Cockburn, Shu and they co-workers in [START_REF] Cockburn | Devising discontinuous Galerkin methods for non-linear hyperbolic conservation laws[END_REF]Karniadakis, Shu, & Cockburn, 2000) DGms distinguish themselves by the definition of fluxes. Among them, there is the method introduced by Dumbser and Käser in [START_REF] Dumbser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes. ii. the three-dimensional isotropic case[END_REF] which employs an Arbitrary high-order DERivation (ADER) scheme which is based on upwind fluxes. The problem with this approach called DG-ADER is that it is dissipative. But as shown in [START_REF] Delcourte | A high-order discontinuous Galerkin method for the seismic wave propagation[END_REF], it is possible to have a conservative scheme by using centered fluxes and employing a Leap-Frog scheme for time integration.

• • • • • • • • • Figure 3.2: h-adaptivity on a 2D mesh in P 1 • • • • • • • • • • • • • • • • • Figure 3.3: p-adaptivity on a 2D
Other schemes have been proposed for second-order wave problems like Interior Penalty DGms ( [START_REF] Grote | Discontinuous Galerkin finite element method for the wave equation[END_REF] for example) which involve a penalty term which makes the scheme stable but it is not always obvious to find the best value of the parameter. More recently, Hybridizable DGms have been introduced, following the pioneering work in [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF]. Regarding the simulation of seismic wave, time-dependent HDGms are fully implicit in time like in [START_REF] Nguyen | High-order implicit hybridizable discontinuous Galerkin methods for acoustics and elastodynamics[END_REF] which can surely explain that HDGms have been developed first for time-harmonic seismic waves. However, we can remark that very recently, some achievements have been carried out for explicit HDG schemes like in [START_REF] Schoeder | Arbitrary high-order explicit hybridizable discontinuous Galerkin methods for the acoustic wave equation[END_REF] for acoustic waves or in [START_REF] Nehmetallah | An explicit hybridizable discontinuous Galerkin method for the 3d time-domain Maxwell equations[END_REF] for elastic waves. It is worth noting that HDGms have been introduced to reduce the number of globally coupled unknowns that are larger with DGms than with Continuous Galerkin approximations. They are DG methods which can be interpreted in terms of Weak Galerkin methods (see the very recent paper [START_REF] Cockburn | The Weak Galerkin methods are rewritings of the Hybridizable Discontinuous Galerkin methods[END_REF]) with the particular feature of involving a stabilization function that ensures the stability and the well-posedness of the formulation. As a conclusion, HDGms are interesting contributions to make the computations faster by reducing the size of the linear systems as compared to DGms but here, since we address the solution of time dependent problems as a preliminary development for inverse problems solutions, we focus on DGms which authorize easily the use of explicit time scheme.

Implementation of DG method for elastodynamics

Here, we address the DG approximation of the elastodynamic problem (2.13) that is meant to be applied in the domain Ω s . We begin with reminding the first-order formulation that we consider for modeling elastic waves into the Earth. Its unknowns are the wavespeed v s and the strain tensor σ which are governed by the following equations: (3.15) where S denotes the source field. We complete the system with initial conditions v s (x, 0) = 0, σ(x, 0) = 0 along with the boundary condition σn s = 0 on the free-surface boundary of the computational domain ∂Ω s . The vector n s denotes the unitary normal vector outwardly directed to Ω s . For constructing the DG formulation of the above problem, we introduce T s,h that defines a collection of Ω s denoted by Ω s,h . T s,h is a set of nonconforming cells K which are triangles or tetrahedra depending on the dimension (2D or 3D). We define Γ s,h the set of all the boundaries of Ω s,h that consists of the external boundary ∂Ω s,h and in addition of all the boundaries of internal cells. In the following, we denote by Γ s,h,out the set of the external boundaries and by Γ s,h,int the set of internal boundaries. Let V s (K) and Σ(K) be two finite subspaces of H 1 (K), typically they are spaces of polynomials. We define a discontinuous approximation (v s,h , σ h ) in V s (K)×Σ(K) as the solution to

       ρ∂ t v s = ∇ • σ, ∂ t σ = C(ǫ(v s )) + S.
           K ρ∂ t v s,h • w s,h = K (∇ • σ h ) • w s,h , K ∂ t σ h : ξ h = K C(ǫ(v s,h )) : ξ. (3.16) for any (w s,h , ξ h ) in V s (K) × Σ(K). Then, using the fact that C(ǫ(v s,h )) : ξ h = Cξ h :
∇v s,h and applying Green formula, we get:

           K ρ∂ t v s,h • w s,h = - K σ h : ∇w s,h + ∂K (σ h n K ) • w s,h , K ∂ t σ h : ξ h = - K (∇ • (Cξ h )) • v s,h + ∂K ((Cξ h )n K ) • v s,h .
where n K denotes the unitary normal vector outwardly directed to K. We move on introducing the discontinuous finite element spaces:

V s,h = {v s ∈ L 2 (Ω s,h ) d | (v s ) |K ∈ V s (K) d ∀K ∈ T s,h }, Σ h = {σ ∈ L 2 (Ω s,h ) d 2 | (σ) |K ∈ Σ(K) d 2 and σ ij = σ ji , ∀K ∈ T s,h }.
The global variational formulation is obtained by summing up each of the contributions of each element K. Then, we have to take new terms into account which correspond to the communications between elements that occur at the internal boundaries Γ s,h,int . It is then necessary to introduce the following quantities in the purpose of defining the numerical fluxes used for gluing the local terms of the variational formulation. Let K 1 and K 2 be two adjacent internal elements. Let u be defined by u K 1 in K 1 and u K 2 in K 2 . We first define [[u]] the jump of u, whose formula depends on the nature of u:

• If u is a scalar we have: [[u]] = u K 1 n K 1 + u K 2 n K 2 ,
which transforms a scalar into a vector

• If u is a vector, we have two possible definitions:

[u] = u K 1 • n K 1 + u K 2 • n K 2 , which
converts a vector into a scalar (such a definition is generally used in the case of acoustic equations) and

[[u]] = u K 1 ⊗ n K 1 + u K 2 ⊗ n K 2
, which transforms a vector into a tensor. We have adopted the notation ⊗ for:

u ⊗ v = u x v x u x v z u z v x u z v z (3.17) in 2D and u ⊗ v =    u x v x u x v y u x v z u y v x u y v y u y v z u z v x u z v y u z v z    (3.18)
in 3D. It is generally used in the case of elastic wave equations.

• If u is a tensor : [[u]] = u K 1 n K 1 + u K 2 n K 2
, which changes a tensor into a vector.

In addition, we define the internal mean value {{u}} which is defined by: {{u}} = 1 2 (u K 1 + u K 2 ) for any internal elements K 1 and K 2 . On the external boundary, we simply have : {{u}} = u while for the jump, we have:

• If u is scalar: [[u]] = un K , • If u is a vector: [u] = u • n K or [[u]] = u ⊗ n K , • If u is a tensor: [[u]] = un K ,
where n K is the unitary normal vector outwardly directed to K and K is an element sharing an edge with the external boundary.

Then by summing up each of the contributions on each element K, we obtain the intermediate DG variational formulation:

             Ω s,h ρ∂ t v s • w s = - Ω s,h σ : ∇w s + Γ s,h [[σw s ]], Ω s,h ∂ t σ : ξ = - Ω s,h (∇ • (Cξ)) • v s +
The remaining terms on the external boundaries can be rewritten as follow:

             Γ s,h,out [[σw s ]] = Γ s,h,out [[σ]] • {{w s }}, Γ s,h,out [[(Cξ)v s ]] = Γ s,h,out [[(Cξ)]] • v s ,
and according to the convention regarding both jump and average associated with edges of the external boundary, we get:

             Γ s,h,out [[σw s ]] = Γ s,h,out (σn s ) • w s = 0 Γ s,h,out [[(Cξ)v s ]] = Γ s,h,out ((Cξ)n s ) • v s .
Then, in order to separate the approximate solution from the test functions, it is convenient to rewrite the fluxes by using the following relations:

       [[σw s ]] = [[σ]] • {{w s }} + {{σ}} : [[w s ]], [[(Cξ)v s ]] = [[Cξ]] • {{v s }} + {{Cξ}} : [[v s ]].
It is worth noting that the solution (v s , σ) to (2.13) is actually continuous inside the computational domain which means that we have:

[[v s ]] = 0, [[σ]] = 0.
Then, we get simplified expressions:

       [[σw s ]] = {{σ}} : [[w s ]], [[(Cξ)v s ]] = [[Cξ]] • {{v s }}.
By plugging these last relations into (3.19), we end up with the full DG formulation:

             Ω s,h ρ∂ t v s • w s = - Ω s,h σ : ∇w s + Γ s,h,int {{σ}} : [[w s ]], Ω s,h ∂ t σ : ξ = - Ω s,h (∇ • (Cξ)) • v s + Γ s,h,out ((Cξ)n s ) • v s + Γ s,h,int [[Cξ]] • {{v s }}. (3.20)
This formulation has been implemented in [START_REF] Boillot | Contributions à la modélisation mathématique et à l'algorithmique parallèle pour l'optimisation d'un propagateur d'ondes élastiques en milieu anisotrope[END_REF]. It turns out that it requires some stabilization terms which act on the formulation as penalties. Basically, the penalized variational formulation we use reads as:

                                         Ω s,h ρ∂ t v s • w s = - Ω s,h σ : ∇w s + Γ s,h,int {{σ}} : [[w s ]] + Γ s,h,int A 1 [[σ]] : [[w s ]] + Γ s,h,int A 2 [[v s ]] : [[w s ]], Ω s,h ∂ t σ : ξ = - Ω s,h (∇ • (Cξ)) • v s + Γ s,h,out ((Cξ)n s ) • v s + Γ s,h,int [[Cξ]] • {{v s }} + Γ s,h,int A 3 [[v s ]][[Cξ]] + Γ s,h,int A 4 [[σ]] • [[Cξ]], (3.21) Each term A 1 , A 2 , A 3 , A 4
, corresponds to a penalty term and there exist different approaches to define them. The most popular formulations are due to Ainsworth et al [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation[END_REF] and [START_REF] Wilcox | A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media[END_REF]. In this work, we use the operators associated with the Ainsworth formulation which has been developed for the acoustic wave equations. It is worth noting that we have extended this formulation to the case of the elastic wave equation. The operators corresponding to the Ainsworth formulation are given in 2D by:

                                 A 1 ([[σ]]) = [[σ]] ⊗ (α) A 2 ([[v s ]]) = α 1 ([[v s ]] xx + [[v s ]] zz ) α 2 ([[v s ]] zx -[[v s ]] xz ) α 2 ([[v s ]] xz -[[v s ]] zx ) α 1 ([[v s ]] xx + [[v s ]] zz ) A 3 ([[v s ]]) = -[[v s ]] : [[α]], A 4 ([[σ]]) = α 3 [[σ]], (3.22)
where

[[v s ]] ij i = x, z, j = x, z are the components of the tensor [[v s ]]. We set α = α • n K
where n K is chosen arbitrarily outwardly or inwardly directed to K. The parameters α, α 1 , α 2 and α 3 are real positive. To the best of our knowledge, there is not any rule for choosing these parameters. Their calibration depends on an empirical process which in turn depends on the simulation parameters.

In the case where α 1 = α 2 = α 3 = 0, the fluxes are said to be centered. The scheme with centered fluxes has the property to preserve the energy and is really adapted for Leap-Frog time-discretization. However, if one of these parameters is non-zero, we obtained upwind fluxes, which implies that the scheme is dissipative and prevents from the use of Leap-Frog discretization. This is why, when it comes to upwind fluxes, we tend to use Runge-Kutta time discretization.

Semi-discrete system

The semi-discrete system associated to the full DG formulation is obtained by choosing the finite dimensional spaces V s (K) and Σ(K). We use to choose V s (K) = P m (K) d and Σ(K) = P m (K) d 2 where P m (K) denotes the space of Lagrange polynomials of degree m defined on the triangle (or tetrahedron) K. However, we can have

V s (K) = Q m (K) d and Σ(K) = Q m (K) d 2
when K is a quadrangle (or hexahedron). An important feature of DG formulation is that we can even have different orders m of approximations and also a combination between P m and Q m .

The basis functions are defined locally, their construction being carried out element per element. Let (ψ K i ) i=1,..,Nm be a basis of the subspace P m (K), where N m is the dimension of P m (K). Here, the dimension N m is m+1 in 1D, (m+1)(m+2)/2 in 2D and (m+1)(m+ 2)(m + 3)/6 in 3D. The set of functions (ψ K i ) i=1..Nm,K∈T s,h allows to construct a basis of V s,h and Σ h by means of cartesian products involving the canonical basis (e j ) j=x,y,z . We then seek the approximate solution in the form:

v s,h|K = j=x,y,z Nm i=1 v K i j (t)ψ i (x)e j .
and

σ h|K = j=x,y,z k=x,y,z Nm i=1 σ K jk i (t)ψ i (x)e j ⊗ e k .
Then, if we choose w s = ψ j e j and ξ = ψ j e j ⊗ e k we can construct the matrices representing the discrete formulation and we get the following matrix system:

                     K∈T s,h   ρ K M K v s,h ∂ t v K s,h + R K v s,h σ K h + L∈T s,h ,K∩L =∅ R K,L σ h σ L h   = 0, K∈T s,h   M K σ h ∂ t σ K h + R K σ h v K s,h + L∈T s,h ,K∩L =∅ R K,L v s,h v L s,h   = 0. (3.23)
where M K ⋆ denotes the mass matrices and R K ⋆ stands for the stiffness matrices, with ⋆ = v s,h or σ h . The important property of these matrices is that they are block diagonal. The matrix

M K v s,h is of size dN m × dN m , with d blocks of size N m × N m . The matrix M K σ h is a d(d + 1)N m /2 × d(d + 1)N m /2 matrix composed of d(d + 1)/2 blocks of size N m × N m .
It is very interesting to observe that the two mass matrices are represented by the entries:

M i,j = K ψ K i ψ K j ∀i, j = 1..N m . (3.24)
which means that there is only one matrix to invert and that the inversion is not necessary for each element because we can go back to the reference element K which is the triangle (0, 0); (0, 1); (1, 0) (Figure 3.4 presents the passage from K to a general element K using the linear transformation F m K ).

(0, 0) (1, 0) (0, 1) F m K (x 3 , z 3 ) (x 2 , z 2 ) (x 1 , z 1 ) Figure 3.4: Passage to the reference element K in 2D
On the reference element K, the entries of the mass matrix become:

M i,j = |detJ F m K | K ψ K i ψ K j = K ψ K i ψ K j ∀i, j = 1..N m .
where

J F m K is the Jacobian of F m K .
On the other hand, the stiffness matrices

R K v s,h and R K σ h are not block diagonal. The matrix R K σ h is a matrix of order d(d + 1)N m /2 × d/N m .
It is dense as the product of the elasticity tensor with the blocks R K i , i = x, y, z that are defined by:

R ipq = K ψ K q ∂ψ K p ∂i . Each block is of size N m × N m . The matrix R K v s,h is of size dN m × d(d + 1
)N m /2 and comprised of the R i blocks defined as follows:

       R K x 0 0 0 R K z R K y 0 R K y 0 R K z 0 R K x 0 0 R K z R K y R K x 0       
Those matrices are also defined using the reference element K and their general term reads:

R ipq = K ψ K q ∂ψ K p ∂i = |detJ F K | K ψp J -T F K ∇ ψq • e i .
The matrices R K,L σ h and R K,L v s,h are the matrices which represent the flux terms between two elements.

R K,L σ h is a matrix of order dN m × d(d + 1)N m /2 and R K,L v s,h is of order d(d + 1)N m /2 × dN m .
They have the same structure than the previous ones, composed of blocks of size N m × N m that read:

               R K,L σ h ,ijkpq = - 1 2 ∂K∩∂L\Γ s,h,out ψ K q e i ⊗ e j n K ψ L p e k , ∀p, q = 1..N m , i = x, y, z j = x, y, z k = x, y, z R K,L v s,h = C(R K,L σ h ) T
As for the mass and stiffness matrices, these matrices can also be written on the reference element.

We then end up with a global discrete system that is written as follows:

M v s,h ∂ t v s,h + R v s,h σ h = 0, M σ h ∂ t σ h + R σ h v s,h = 0. (3.25)
with the following details:

(M v s,h ) K,L = ρ K M K v s,h if L = K, 0 if L = K, (3.26) (R σ h ) K,L =        R K σ h + R K,K σ h if L = K, R K,L σ h if L = K and ∂L ∪ ∂K = ∅, 0 else, (3.27) (M σ h ) K,L = ρ K M K σ h if L = K, 0 if L = K, (3.28) (R v s,h ) K,L =        R K v s,h + R K,K v s,h if L = K, R K,L v s,h if L = K and ∂L ∪ ∂K = ∅, 0 else.
(3.29)

Validation of DGm on structured quadrangle meshes

In the purpose of motivating the coupling of DG finite elements with spectral elements, we have implemented the previous DG formulation on quadrangle meshes, considering that spectral elements are particularly adapted to quadrangle meshes. In what follows, we provide a validation of our implementation which is a new feature for Elasticus. It is worth noting that we have addressed both the acoustic and the elastic systems. For more information on the acoustic system we consider, we refer the reader to the section related to the spectral element method.

The simulation domain Ω s is a 3000 m×3000 m square (see Figure 3.5 for the general settings in particular the position of the source) and we consider two different configurations:

• A homogeneous acoustic domain with P -wavespeed equal to 1000m.s -1 and ρ equal to 1kg.m -3 .

• A homogeneous elastic domain with V P = 1000m.s -1 and V S = 500m.s -1 and ρ = 1kg.m -3 . The validation consists in comparing the numerical solution obtained with Elasticus with an analytical one computed with the software Gar6more developed in the team Magique3D (gar6more2d.gforge.inria.fr). To assess possible differences between both computations, we picture seismograms which represent the time variation of the velocity in the x direction in a given point of the computational domain. For the seismograms, the solution for both Gar6more and Elasticus is recorded thanks to a receiver placed inside the computational domain. Here, the receiver is placed at point (1000,2000). The square is covered with a collection of quadrangles and we consider two grids of 10000 cells and 22500 cells. The simulations are performed with a space approximation of order 3 and the time integration uses a second order Leap-Frog scheme whose CFL number is accurately computed with the power iteration. All computations are performed in parallel using OpenMP with eight threads. In each case, the wave propagation is generated by a P -Ricker point source located at the center of the domain with f peak = 10Hz and t peak = 0.12s and we apply a free-surface boundary condition on each side of the propagation domain. The DG scheme includes Ainsworth penalization term with coefficient α equal to 0.5 and the other parameters set to zero. We illustrate the simulations with two kinds of results: some snapshots which give an overall behavior of the waves in each case and seismograms obtained with penalized DGm and superimposed with the analytical solution computed with Gar6more. Figures 3.8 and 3.9 picture the seismograms in the acoustic case for the two different meshes. We observe that the two curves are very close and this qualitative evaluation is confirmed by computing the relative L 2 -error displayed at formula (3.30) which is 1.3e -3 (0.13%) with a mesh of 10000 cells and 5.8e -4 with a mesh of 22500 cells. To compute this L 2 -error, we take the numerical solution noted u Elasticus and the analytical one noted u Gar6more and we compute the following quantity:

Err L 2 = n i=1 (u Elasticus (x 0 , i) -u Gar6more (x 0 , i)) 2 n i=1 u Gar6more (x 0 , i) 2 , (3.30)
where n is the number of time steps and x 0 the position of the receiver we consider.

To get more insight on the accuracy of the DG simulations, we have pictured Figures 3.6 and 3.7 first and Figures 3.10 and 3.11 next. All the figures show that the propagation of waves is correctly reported, both in the acoustic and the elastic case. In the elastic domain we can observe the absence of S-waves. This is due to the P -wave source we generate and the homogeneous nature of the middle. However, the S-waves appear when the wave is reflected on the boundary as seen in Figures 3.10 and 3.11. As a conclusion,we have introduced the DGm that is used in this work to solve wave problems formulated as first-order systems and we have presented some numerical results of DGm implemented with quadrangle meshes. This is a feature of Elasticus that did not exist at the beginning of this study. It has been developed to provide a way to compare DGm with spectral finite elements that are more adapted to quadrangles than triangles. The idea of comparing these two approaches comes from our commitment to couple DGm with SEm. The latter is the purpose of the next section.

Mixed-primal formulation of DG

We have introduced the DG formulation that we use to apply for solving wave problems and it is done by many other authors [START_REF] Delcourte | A high-order discontinuous Galerkin method for the seismic wave propagation[END_REF][START_REF] Ainsworth | Dispersive and dissipative behaviour of high order discontinuous Galerkin finite element methods[END_REF] because it has the advantage of being symmetrical. However, we aim at coupling DGm with SEm and it turns out that the coupling is not that obvious when considering this formulation. This is why we introduce here another formulation, which is strictly equivalent to the former formulation (3.20).

             Ω s,h ρ∂ t v s,h • w s,h = - Ω s,h σ h : ∇w s,h + Γ s,h,int {{σ h }} : [[w s,h ]], Ω s,h ∂ t σ h : ξ h = - Ω s,h (∇ • (Cξ h )) • v s,h + Γ s,h,out ((Cξ h )n s ) • v s,h + Γ s,h,int [[Cξ h ]] • {{v s,h }}.
(3.31)

We consider now the second equation to which we apply an integration by parts:

Ω s,h ∂ t σ h : ξ h = -- Ω s,h (Cξ h ) : ∇v s,h + Γ s,h,out ((Cξ h )n s ) • v s,h + Γ s,h,int [[Cξ h v s,h ]] + Γ s,h,out ((Cξ h )n s ) : v s,h + Γ s,h,int [[Cξ h ]] • {{v s,h }}.
Then, the terms related to the external boundary cancel. Moreover, we have the following relation:

[[Cξ h v s,h ]] = [[Cξ h ]] • {{v s,h }} + {{Cξ h }} : [[v s,h ]].
which leads to the mixed-primal DG formulation: which is equivalent to the initial one for which we dispose of stability and convergence results [START_REF] Delcourte | A high-order discontinuous Galerkin method for the seismic wave propagation[END_REF]. We will see later on that this formulation shares all what is needed with the SE formulation to obtain a natural coupling of these two methods.

             Ω s,h ρ∂ t v s,h • w s,h = - Ω s,h σ h : ∇w s,h + Γ s,h,int {{σ h }} : [[w s,h ]], Ω s,h ∂ t σ h : ξ h = Ω s,h (Cξ h ) : ∇v s,h - Γ s,h,int {{Cξ h }} : [[v s,h ]].

Spectral element formulation

In this section, we introduce the spectral element method that has been widely used by geophysicists for solving problems of seismology. Prior to any numerical development leading to performance assessent, it distinguishes itself from DGm at the mesh level. Indeed, its development is more suitable for quadrangles/hexahedra than triangles/tetrahedra while DGm can be applied to any type of meshes including those with non conformal elements.

Overall setting of SEm

The Spectral Element (SE) method has been introduced by Patera and Maday for fluid dynamics in [START_REF] Patera | A spectral element method for fluid dynamics: laminar flow in a channel expansion[END_REF][START_REF] Maday | Spectral element methods for the incompressible navier-stokes equations[END_REF]. Regarding seismic applications, it has been used first by Priolo, Carcione and Seriani in 1994 [START_REF] Priolo | Numerical simulation of interface waves by high-order spectral modeling techniques[END_REF]. At the beginning, it employed Chebyshev polynomials [START_REF] Faccioli | Spectral-domain decomposition method for the solution of acoustic and elastic wave propagation[END_REF] and soon after, Komatitsch and Vilotte proposed to use Lagrange polynomials (Komatitsch & Vilotte, 1998). It is now widely used by Geophysicists since the tremendous work of Komatitsch and his co-authors [START_REF] Komatitsch | Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2 D et 3 D en milieu hétérogène[END_REF]Komatitsch et al., 1998;[START_REF] Komatitsch | Wave propagation near a fluid-solid interface: A spectral-element approach[END_REF]Komatitsch et al., 2000a;[START_REF] Komatitsch | Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles[END_REF] who have largely contributed to both the development and maintenance of the software package Specfem (https://geodynamics.org/cig/software/specfem3d/). Specfem has demonstrated a great ability in accurately and fastly reproducing wave propagation in the framework of 3D global and local seismology. Today, besides Specfem, Salvus (https://www.pasc-ch.org/ projects/2017-2020/salvus/) is a very efficient open-source platform under development that perform simulation of large-scale waveform modelling and inversion which employs SEm as well and proposes frequency-domain solvers while Specfem focuses on time-domain problems.

SEm couples the flexibility of standard Finite Elements and the accuracy of pseudospectral methods. By using Gauss-Lobatto-Legendre points as degrees of freedom and applying Gauss-Lobatto-Legendre quadrature method for the computation of the Lagrange polynomials integrals involved by the variational formulation of the boundary value problem, SEm distinguishes itself from other continuous finite element approximations by having a diagonal mass matrix. This is a clear advance as compared with standard continuous approximations by allowing explicit time integration without involving any mass lumping process which is well-known for causing a decrease of the order of convergence.

• • • • • • • • • • • • • • • • • • • • • • • • • Figure 3.14: Position of the Gauss-Lobatto points in Q 4
SEms use tensor-product Lagrange interpolants within each element and the nodes of these shape functions correspond to Gauss-Lobatto-Legendre points which are the zeros of Legendre polynomials. For any smooth function, it has been proven that its interpolant converges exponentially fast as the order of the interpolant is increased, which justifies the method is qualified as spectral. Moreover, this simplifies the resolution of the matrix system because the inversion of the mass matrix is direct and thus at no cost. The quadrature also makes the stiffness matrix sparse and thus contributes to speed up the computations.

In terms of wave propagation, SEm had been introduced for the elastic and acoustic 2D case by [START_REF] Seriani | High-order spectral element method for acoustic wave modeling[END_REF][START_REF] Seriani | High-order spectral element method for elastic wave modeling[END_REF]) using a Chebyshev polynomial basis. In their studies, they conclude that a low number of grid points per wavelength is enough to ensure stability even at high order. Next, they extended their study to the simulation of acoustic wave simulation in heterogeneous media (Seriani & Priolo, 1994;[START_REF] Seriani | Modelling waves in anisotropic media by a spectral element method[END_REF], where they used a Fourier analysis to confirm the robustness of SEm regarding numerical dispersion.

The elastodynamics system in 2D and 3D heterogeneous media was addressed in the thesis of Komatitsch, [START_REF] Komatitsch | Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2 D et 3 D en milieu hétérogène[END_REF]. This work has been followed by a large number of publications. In particular, the efficiency of SEm for applications dealing with seismic problems in complex domains was investigated in (Komatitsch & Vilotte, 1998;Komatitsch et al., 1998) for the elastic case. The extension to the elasto-acoustic case was done in [START_REF] Komatitsch | Wave propagation near a fluid-solid interface: A spectral-element approach[END_REF] where the SEm is introduced to deal with the interface between an acoustic medium and an elastic medium. For anisotropic media, we can quote the investigation work of Komatitsch et al. in (Komatitsch et al., 2000a) where difficulties of some numerical methods, like the Finite Differences, are displayed to convince that Lagrange interpolants are really necessary in anisotropic cases like VTI. Indeed, with the use of other bases like Chebyshev polynomials, the mass matrix is no longer diagonal. One of the main disadvantage of the classical SEm is that the mesh is composed of quadrangles in 2D, hexehadra in 3D, which have difficulties to respect the features of a geological medium such as velocity contrasts or more generally heterogeneities. Triangular or tetrahedral meshes are much more flexible and one approach consists in first coupling SEm with standard finite elements to gain more flexibility in some areas of the domain. Such an idea was formerly explored in [START_REF] Moczo | Hybrid modeling of p-sv seismic motion at inhomogeneous viscoelastic topographic structures[END_REF] to increase the flexibility of finite differences by a coupling with standard finite elements. For instance, the mortar method [START_REF] Lahaye | Hybrid finite element-spectral element approximation of wave propagation problems[END_REF][START_REF] Casadei | Implementation of a 3d coupled spectral element solver for wave propagation and soil-structure interaction simulations[END_REF] provides an interesting methodology for nonconforming matching between spectral element and finite element meshes but it turns out that it increases significantly the computational costs. Some works have been done to develop SEm on triangle and tetrahedral grids. The first idea that comes to mind is the construction of Gauss Lobatto-Legendre points for triangles. Unfortunately, except in very particular situations like segment [-1, 1], it is unclear that such points exist, and if they do, their construction for numerical use should be difficult. Few results on that approach for the triangle are summarized in [START_REF] Lyness | A survey of numerical cubature over triangles[END_REF]. The work of Karniadakis and Sherwin in [START_REF] Sherwin | A triangular spectral element method; applications to the incompressible Navier-Stokes equations[END_REF] is particularly interesting. In this work, they use the orthogonal basis and the polynomials defined by Dubiner in [START_REF] Dubiner | Spectral methods on triangles and other domains[END_REF] where the inner product is replaced by a warped product, designed for accurately approximating the integrals. Unfortunately, in this kind of grids the warped product is oversampled: it requires twice as many points as there are degrees of freedom and as a consequence, the mass matrix is no longer diagonal. To avoid this issue, Taylor and Wingate have proposed to use a new set of points in [START_REF] Taylor | The Fekete collocation points for triangular spectral elements[END_REF], i.e. the so-called Fekete points. The Fekete points are a natural generalization of the points used by the SEm on quadrangles. They enable a conforming matching between triangles and quadrangles and the important properties of SEms are kept including the diagonal mass matrix. Such an interesting property is satisfied because Fekete points located at the edges of triangles (faces of tetrahedra) coincide with Gauss-Lobatto-Legendre points placed at the edges of quadrangles (faces of hexahedra). By this way, mortar elements are useless for the coupling, triangles and quadrangles are conform and the resulting mass matrix is diagonal. In [START_REF] Komatitsch | Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles[END_REF], the use of Fekete points is applied with success to the modeling of wave propagation in 2D elastic media. They employ Dubiner polynomials that are coupled with Fekete points in order to keep the global mass matrix diagonal in the part of the mesh composed of triangles. Numerical experiments show that triangles are slightly less accurate than quadrangles. In particular, the number of points per wavelengths must be increased to keep the same level of accuracy with triangles and it has consequence on the time step by respecting the stability condition of Courant-Friedrich-Levy. Moreover, the use of tensor-product is no longer possible on triangles which results in increasing the number of operations for computing a derivative. In particular, the authors show that the computation in the triangular area of the mesh is more expensive than in quadrangle one with a ratio of ≃ N/2 where N is the degree of the polynomial basis and the extra cost should be higher in 3D cases. However, they observe that the computational extra cost is compensated by a gain of flexibility provided by the use of two types of elements.

The conclusions of [START_REF] Komatitsch | Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles[END_REF] indicate that there is a clear interest in bringing some more flexibility to SEm and that can be done by working with hybrid meshes composed of quadrangles and triangles. The efficiency of SEm on quadrangles and hexahedra is undeniable but on triangles, the results on [START_REF] Komatitsch | Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles[END_REF] indicate some issues as computational extra costs and decreasing of the time-step value. This is then a clear motivation to develop another coupling by adopting another point of view. Rather that seeking to bring more flexibility to SEm, we propose to couple SEm on quadrangles (or hexahedra) with a method which has demonstrated its efficiency on triangles (or tetrahedra).

Implementation of the Spectral Element method

At the beginning of this work, the team Magique-3D developed and maintained advanced software package mainly based on DGms. Besides all the advantages provided by DG approximations, the team was particularly interested in using triangular or tetrahedral meshes to ensure a high level of flexibility for the approximation of the geological layers. Nevertheless, we have many examples where the use of conforming grids is not necessary, for instance in homogeneous media and the excellent performances of Spectral Element method (SEm) for solving wave problems offer a great opportunity to improve the performance of our software in the purpose of a coupling. That is why we have decided to work on the coupling between the DG-based software and a SE-based code that has been developed first in the framework of this thesis. We could have used Specfem, which has demonstrated great capacities for solving wave problems but the team has an agreement with the industrial group Total which develops and maintains its own computational platform. Given this, it turned out more relevant to develop our own SE code respecting the industrial environment of development, in particular when deciding data structures compatible with further coupling. This section deals with this part of development which has been done from scratch. To introduce the SEm, we will focus on the acoustic system which offers the simplest configuration, one of the two unknowns being scalar. But it is worth noting that the method has also been implemented for elastic wave problems with a very similar approach, which does not justify to be presented in details as well. To illustrate this part of achievement, we will simply provide some numerical results in the elastic case.

SEm variational formulation

Let Ω f be the fluid domain in which we solve the acoustic wave equation. We consider the first order formulation of the acoustic wave equation having the velocity vector v f and the pressure p as unknowns:

                               ρ∂ t v f = -∇p in Ω f × (0, T ), 1 c 2 0 ρ ∂ t p = -∇ • v f + S in Ω f × (0, T ), p(x, 0) = 0, v f (x, 0) = 0 in Ω f , v f • n f = 0 on ∂Ω f × (0, T ).
(3.32)

Here T denotes the final time of computation, ρ denotes the density of the fluid and c 0 stands for the velocity of the waves in the fluid. ∂Ω f is the boundary of the fluid domain.

It is supposed regular enough to define a unitary normal vector n f which is outwardly directed to Ω f . For simplicity, we assume that the boundary is polygonal or polyhedral, which implies that we can define a normal vector locally on each edge of the boundary. S is the source term which depends both on time and space. We seek a solution (v f , p) in H(div, Ω f ) × H 1 (Ω f ). For that purpose, we consider that for any pair of sufficiently smooth test-functions (w f , q), we have:

             Ω f ρ∂ t v f • w f = - Ω f ∇p • w f , 1 c 2 0 ρ Ω f ∂ t pq = - Ω f ∇ • v f q + Ω f
Sq.

(3.33)

Then, we apply a Green formula on the second equation of (3.33):

1 c 2 0 ρ Ω f ∂ t pq = Ω f ∇q • v f + < v f • n f , q > -1/2,1/2 + Ω f Sq (3.34)
where <, > -1/2,1/2 denotes the duality brackets of trace spaces H -1/2 (∂Ω) for H(div, Ω) and H 1/2 (∂Ω) for H 1 (Ω). The boundary term vanishes according to the boundary condition on the velocity and we end up with the following variational formulation:

Find (v f , p) ∈ (L 2 (Ω f ), H 1 (Ω f )) such that, for any couple (w f , q) ∈ (L 2 (Ω f ), H 1 (Ω f )), we have:        Ω f ρ∂ t v f • w f = - Ω f ∇p • w f , 1 c 2 0 ρ Ω f ∂ t pq = Ω f ∇q • v f + Ω f Sq (3.35)

Semi-discretization

To construct the semi-discrete system, we first introduce T f,h , which defines a paving of Ω f denoted by Ω f,h . T f,h is a set of quadrangle/hexahedral cells and K stands for an element of T f,h . The set of element edges common with ∂Ω f is denoted by Γ f,h . Then we define the finite element spaces:

V f,h = {v f ∈ (L 2 (Ω f,h )) d | (v f ) |K ∈ (Q l (K)) d , ∀K ∈ T f,h }, P f,h = {p ∈ H 1 (Ω f,h ) | (p) |K ∈ Q l (K), ∀K ∈ T f,h }.
In the above definition, Q l (K) is the space of polynomials with degree inferior or equal to l. Then, following [START_REF] Diaz | Approches analytiques et numériques de problèmes de transmission en propagation d'ondes en régime transitoire. application au couplage fluide-structure et aux méthodes de couches parfaitement adaptées[END_REF], we seek for an approximate solution

(v f,h , p h ) ∈ V f,h × P f,h satisfying:              Ω f,h ρ∂ t v f,h • w h = - Ω f,h ∇p h • w h , Ω f,h ∂ t p h q h = Ω f,h ∇q h • v f,h + Ω f,h Sq h , (3.36) for any (w h , q h ) in V f,h × P f,h .
Let K be the reference element (unit square in 2D, unit cuboid in 3D). We introduce a linear transformation F m K which defines a basis function on a given element K m from the basis functions defined in the reference element K (Figure 3.15 depicts the passage from

K to K m ). (0, 0) (1, 0) (0, 1) (1, 1) F m K (x 4 , z 4 ) (x 3 , z 3 ) (x 2 , z 2 ) (x 1 , z 1 )
Figure 3.15: Passage to the reference element K in 2D

Here, we work with the polynomial space Q l , whose construction requires using (l +1) d points. In the case of SEm, we choose the Gauss-Lobatto points associated with the reference element and denoted by ηi , i = 1, .., (l + 1) d . Then, the basis functions defined in K are Lagrange polynomials constructed by solving the equations φj (η i ) = δ i,j , (3.37)

where δ i,j stands for the Kronecker symbol, i, j = 1, .., (l + 1) d . From this, we have the basis functions of the element K m thanks to the expression

ϕ m j = φj • F m K , j = 1, .., (l + 1) d . (3.38)
Then, the approximate v f,h is sought in the form

v f,h = j=x,y,z ndof i=1 v i j (t)ϕ i (x)e j ,
where e j j=x,y,z stands for the canonical basis.

If the mesh is composed of L nodes N j , we define Φ j , j = 1...L as the Lagrange polynomials constructed from the nodes. Functions Φ j satisfy the following properties:

• the support of Φ j is composed of the set of all quadrangles/hexahedra K which contain the j-th degrees of freedom;

• Φ m j = Φj • F m K , j = 1, .., L and Φj ( Ni ) = δ i,j .
Then, the interpolate of the pressure p h is given by:

p h (x, t) = (l+1) d i=1 p i (t)Φ i (x),
where Φ j , with j = 1, .., L, are the polynomial basis functions previously introduced.

It is worth noting that, by construction of the finite element spaces, the approximate velocity is sought in a space of discontinous functions while the approximate pressure is continuous. Figure 3.16 pictures both the local degrees of freedom of an element K (in red) and the global degrees of freedom (in black). The wave velocity v f,h , as a discontinuous function, admits a local representation involving the red points in each element. They correspond to the coefficients used in the decomposition with the basis functions ϕ i that are Lagrange polynomials computed with the Gauss-Lobatto points. Regarding the

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 3
.16: Elements with local Gauss-Lobatto points(in red) and degrees of freedom on the mesh (in black) in Q 3 . Note that the local points inside the elements coincide with the global one.

approximate p h , it is defined globally in the mesh from discrete values computed at the degrees of freedom in black.

By choosing basis functions as test functions, we end up with the matrix system

ρM v f,h ∂ t V f,h + R v f,h P f,h = 0, M p h ∂ t P f,h + R p h V f,h = S h , (3.39)
where

• V f,h is the discrete wavespeed vector whose components are obtained from the decomposition of the v f,h in the basis φ i ,

• P f,h is the discrete pressure vector whose components are obtained from the decomposition of the p h in the basis Φ i ;

• M v f,h and M p h are the mass matrices;

• R v f,h and R p h are the stiffness matrices. By construction, we have

R v f,h = R T p h ;
• S h is the discrete source term decomposed in the basis Φ i .

As previously mentioned, we compute the integral coefficients of each matrix by using Gauss-Lobatto quadrature formulas. For instance, let (η i ) be the Gauss-Lobatto-Legendre points and (ω i ) the corresponding quadrature weights. Then, we have:

(M v f,h ) i,j ≃ e∈supp(ϕ i )∩supp(ϕ j ) (l+1) d k=1 ω k ϕ i (η k )ϕ j (η k ).
By construction, we have that ϕ i (η k ) = δ i,k , where δ i,k stands for the Kronecker symbol. The mass matrix is thus represented by the coefficients:

(M v f,h ) i,j ≃ e∈supp(ϕ i )∩supp(ϕ j ) ω i δ i,j .
The same conclusion holds for the other mass matrix. On the reference element K, the general terms of the mass matrices become:

M ij ≃ e∈supp(ϕ i )∩supp(ϕ j ) ω i δ i,j ≃ (l+1) d k=1 |detJ F K |ω i δ i,j .
where

J F K is the Jacobian of F K.
Regarding the stiffness matrices, the use of Gauss Lobatto formulas for computing their entries simplifies the calculations but does not change the structure of the matrices.

Validation of SEm

As for DGm, we validate here the implementation of SEm on structured quadrangle meshes. The validation is achieved by considering the same test cases than for DGm validation, which means that we have exactly the same data (point source, physical parameters) and use the same time scheme : the Leap-Frog scheme. Figures 3.17 (1000,2000) shows that the numerical solution follows perfectly the profile of the analytical solution. We have computed the relative L 2 -error between the two solutions: it is equal to 8e -3 for the grid with 10000 cells and 3.6e -3 for the one with 22500 cells.

We reproduce the same experiment in the case of an elastic domain. Figures 3.21 and 3.22,completed with Figures 3.23 and 3.24, show that the same conclusions hold.

Remark 1: Again in the elastic case, the S-wave only appears when the wave hits the boundary due to P -wave source and the homogeneous nature of the domain.

Conclusion

We have implemented and validated the SEm for solving the acoustic wave equation formulated as a first order system having the velocity and the pressure as unknowns. The approximate pressure is continuous while the velocity is discontinuous. Despite our will to apply DGm on triangle/tetrahedral meshes, we have also implemented DG for quadrangle meshes. By this way, we dispose of a wide variety of numerical schemes which will be used in what follows as proof of concept of our coupling strategy.

Comparison of DGm with SEm

In this section, we compare DGm and SEm when applied in structured quadrangle meshes. The DG scheme involves centered fluxes which ensures the scheme to be conservative. The .24: Seismograms of Gar6more solution (in blue) and numerical solution using SEm (in red) applied in an elastic domain discretized with 22500 quadrangle cells comparison is performed in terms of accuracy and CPU time. We use the Leap frog scheme for time integration and the computations are carried out in parallel using OpenMP on 8 threads. We begin with computations in homogeneous 2D domains which have already been introduced in the previous sections and we consolidate our study by providing an example of 2D heterogeneous domain.

Comparison in terms of CFL number

Explicit time integrations require following a CFL condition which ensures the stability of the numerical process. Its value depends on the space discretization and on the scheme used for the time integration. In practice, the optimal CFL number is computed using the power iteration method algorithm which computes the highest eigenvalue ρ(A) of the matrix A governing the final linear system. Once this eigenvalue is computed, the CFL number is obtained as the square root of ρ(A) -1 multiplied by a constant related to the time scheme. Here, the constant associated to Leap-Frog scheme is equal to 2 and we use to employ 1.99 to strengthen the validity of the CFL condition. The comparison is based upon the L 2 relative error which is computed by using an exact solution obtained with Gar6more software and the CPU time.

Table 3.5 pictures the CFL number and the computational time for the two numerical methods in the case of an acoustic domain paved with 10000 cells previously described in chapter 3. To begin with, the DG scheme does not involve any penalization term, which is the simplest DG formulation of the problem.

CFL

L2-error CPU-time Nb of time steps DG 1.56e-3 0.941 4.98 629 SEM 2.45e-3 7.903e-3 0.80 409 Table 3.5: Acoustic case, DG scheme without penalization, DGm and SEM of order 3, 10000 quadrangle cells.

As expected, the two CFL numbers are not the same but they are not so far from each other. The SEm clearly provides more accuracy than DGm as pictured in the third column which indicates that the DG error is about 94%. The SEm CPU time is really smaller than the one of DGm as it is 6 times lower. This first experiment is thus clearly favouring the SEm. Regarding DG accuracy, it is worth noting that DGm is not penalized and this should explain why the numerical error is so high. The lack of penalization is illustrated in Figure 3.25, which depicts a P -wave clearly polluted by numerical noises. This observation is confirmed globally when computing the relative L 2 -error thanks to the analytical solution computed with Gar6more which is almost 100%. In addition to the lack of penalization, such a high numerical error can also be explained by the use of structured cells. Indeed, computations in structured meshes can be polluted by resonance phenomena that bring out the numerical noise. Hence, in the following, we will systematically consider a DGm involving Ainsworth penalization terms and in all of our experiments, we consider a penalization coefficient α equal to 0.5 while the other penalization coefficients are set to 0. In Table 3.6, we can see that adding a penalization term in the DG formulation restores the accuracy of the DG solution and we can even see that the DG error is almost six times lower than SE one when applied in the coarser mesh while it is more than ten times lower in the finer mesh. However, the computational costs for DG method are increased. This is mainly due to the modification of the CFL number which is smaller as a consequence of introducing a penalization term. Just as was formerly seen in other works [START_REF] Ventimiglia | Schémas numérique d'ordre élevé en temps et en espace pour l'équation des ondes du premier ordre. application à la Reverse Time Migration[END_REF], the addition of penalization terms actually modifies the spectrum of the matrix A such that the value of ρ(A) increases. This contributes then to decrease the CFL number as indicated in We also display the results we have obtained in the elastic case in Table 3.7 and we can see that the same conclusions hold.

As a conclusion, in terms of CFL number and thus CPU time, the SEm is more advantageous than DGm when implemented in structured quadrangle cells. However, we see that when using the optimal CFL number for each scheme, the L 2 -relative error of DGm is always lower than the one of SEm.

Comparison in terms of relative error

We have seen in Table 3.6 that, for a given mesh and with the use of the optimal CFL number, the error of penalized DGm is smaller than the one of SEm. However, the computational cost of DGm is always higher than the one of SEm. We must then keep in mind that the objective of this work is the coupling of both methods and we aim at providing a numerical approach that monopolizes less computational resources while ensuring a good accuracy of the simulations. It would thus be interesting to adjust the numerical parameters in order to get a balance between numerical error and computational time. This idea led us to draw Figure 3.26 which displays the evolution of the relative error as a function of the approximation order. The underlying simulations have been performed by using the optimal CFL number of each method. For instance, we get the information that for the considered example, DGm and SEm deliver a solution with a relative error of about 1.3e -3 when the approximation degree is 3 for the first and 5 for the latter. In Table 3.8, we have then displayed the values corresponding to the DGm of order 3 and the SEm of order 5. We see that the CPU time is still higher for DGm than SEm even if the order of SEm is higher. We must realize here that this experiment is quite lucky in the sense that both methods have the same optimal CFL, as pictured in the second column of Table 3.8. This could indicate that we can do just as well by employing SEm of lower order with the optimal CFL of DGm. We just have to compute the optimal CFL number of lower order SEm to be sure that the one of DGm is correct. For instance we have computed the optimal CFL of the SEm at order 3 that is about 1.06 e-3 and thus higher than the one of DGm. In Table 3.9, we then see that we reach the same accuracy for the SEm at order 3 as for the DGm at order 3 providing the SEm is used with the optimal CFL number of SEm at order 5.

To complete this study, we have compared the L 2 -error of DGm at order 3, SEm at order 3 and SEm at order 5 when using 20, 40, 60, 80 and 100% of the optimal CFL of each method. The results are depicted on Figure 3.27 and confirmed that when the CFL of DG at order 3 is close to that of the SEm at order 3, the error are also close to each other. 2.12 1000 Table 3.9: Penalized DGm of order 3 and SEm of order 3, acoustic domain, 10000 quadrangles

Comparison in terms of meshes

According to the results of the previous subsection, we use here the optimal CFL number computed for the DGm of order 3 and we employ it both for DGm and SEm of order 3.

Table 3.10 depicts the results for the acoustic case in each mesh and the elastic case is displayed in Table 3

Bi-layered acoustic case

In this section, we address the case of a heterogeneous domain composed of two layers separated by an irregular (staircase) flat interface. We aim at confirming that, here again, the SEm is still more advantageous than DGm in an heterogeneous acoustic medium paved with structured quadrangle cells. The domain of study we consider is depicted in Figure 3.28. In terms of physical parameters, the acoustic wavespeed in the top layer is 1000m.s -1 and 2000m.s -1 in the bottom layer. For both layers, the density ρ is 1kg.m -3 .

• src 1980m 1500m Figure 3.28: Bi-layered acoustic domain with a stair-like interface.

This boundary domain is a square and the computational area is paved with 10000 squares of side 3m. The origin of the system of coordinates is located in the bottom left corner. Then, the source term is a second order Ricker point source placed at point (765,2450) computed using the formula given in Chapter 2 on page 28 with a peak frequency equals to 5Hz and a time delay equals to 0.24s. The boundary condition is a free surface condition which is set on the whole boundary. We employ DGm and SEm of order three coupled with the Leap Frog scheme in time and, as in the previous comparisons, we base our analysis on the computational time of the simulation and the value of the relative L 2 -error computed from a reference solution constructed with Elasticus using DGm of order six in space. Indeed, the geometry of the interface prevents the domain of propagation from testing the hypotheses required for using Gar6more software. As in the previous section, the optimal CFL number is computed with the power iteration method. Table 3.12 pictures the results. We observe that the CPU-time of SEm is still very short in comparison with the DG-one with a factor of nearly sixteen between the two values. However, the SEm error is 20 times greater than the DGm error. Hence, in order to balance both errors, we have relaunched the simulation with the SEm for the same CFL number than the one of DGm, which means that the computations are performed with a smaller time step. The corresponding comparison is displayed in Table 3.13. As expected, the use of a smaller CFL number for applying SEm contributes to reduce the error that becomes nearly equal to the one of DGm and obviously to increase the computational time. But the CPU-time of SEm remains four time smaller than the one of DGm. Hence, this example shows that the same conclusions hold in the case of a heterogeneous domain.

CFL

Conclusion

We have compared DGm and SEm applied in structured quadrangle meshes. The analysis is focused on the case of waves propagating in acoustic domains, even if the methods have also been developed for elastic domains. The numerical experiments we have performed show that the SEm is clearly faster than DGm on structured quadrangle cells. At first sight, we have applied the numerical methods with their respective optimal CFL number. In that case, we have observed that the SEm is less accurate than the DGm. For the purpose of balancing the two errors, we propose to use the SEm with the optimal CFL number of DGm. Then, the computational cost of SEm is increased while remaining smaller than the one of DGm and the level of accuracy is maximized. This result opens the possibility of coupling efficiently the two methods. All these comparisons have been performed with Leap-Frog time scheme. It is worth noting that some numerical noise can pollute the simulations with DGm. This numerical pollution can be controlled by adding a penalization term as illustrated in this chapter. However, in [START_REF] Ventimiglia | Schémas numérique d'ordre élevé en temps et en espace pour l'équation des ondes du premier ordre. application à la Reverse Time Migration[END_REF], it has been shown that if the mesh is composed of structured triangle cells, the addition of a penalization term may be ineffective in the sense that the amplitude of the numerical noise is decreased but not enough to guarantee a good level of accuracy. To illustrate the phenomenon, we have performed a simulation in an acoustic squared domain of side 3000 meters. The P -wavespeed is equal to 1000m.s -1 and the density ρ is equal to 1kg.m -3 . This domain is paved with 15000 structured triangles and we apply DGm of order 3 with Ainsworth penalization (α = 0.5 and the other parameters equal to zero). The propagation is initiated with a Ricker point source placed at the center of the domain with a central frequency f peak equals to 20Hz and a time delay equals to 0.04s. Numerical results are depicted in Figure 3.31 where we can see that even with penalization term, the numerical noise actually pollutes significantly the numerical solution.

For this reason, we have decided to move on with Runge-Kutta time-schemes which have proven to be stable when using upwind fluxes. It is worth noting that in this case, the scheme is now dissipative.

Conclusion

We have presented and studied the two numerical methods that form the base of our work by focusing on their application in quadrangle meshes. For the sake of conciseness of this manuscript, we have limited the introduction to the elastodynamics for DGm and acoustics for SEm. However, each numerical method has been implemented in Elasticus both for elastic and acoustic waves in 2D and 3D. This chapter contains some numerical results in 2D which aim at validating the development of DGm and SEm in structured quadrangle cells. We have then developed a numerical comparison of both methods. We have observed that:

• For a given order of approximation in space and the use of a Leap-frog scheme, the CFL number of SEm is always lower than the one of DGm. Hence it could have been interesting to introduce local time stepping in order to apply a time step that is adapted to each numerical method.

• When employing the optimal CFL number for each method, the DGm delivers a solution that is clearly more accurate that the solution obtained with SEm. This difference can be attenuated by employing the SEm with the optimal CFL number of DGm. Then the CPU time of SEm increases but is still slearly lower than the one of DGm.

Given these two previous items, we can conclude that it is more relevant to use DGm and SEm with the same time step providing it is defined by the DGm rules.

Chapter 4

DG and SEm coupling strategy

This chapter deals with the coupling of the SEm and the DGm that have been introduced and analysed in the previous chapter. The idea of coupling different finite element methods has been addressed in many works and each approach distinguishes itself from others by the matching condition involved for the coupling. For instance, the mortar elements introduced in [START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method[END_REF] have demonstrated a high efficiency for coupling finite element methods. The mortar method involves a matching condition on the interface between two different subdomains of the mesh. Then a numerical method is applied, possibly specific to each subdomain. This approach usually leads to a nonconforming method providing optimal approximation. The coupling of SEm with another numerical method has already been addressed in the past. For instance, Bernardi and her collaborators [START_REF] Bernardi | Coupling finite element and spectral methods: First results[END_REF] coupled SEm with a finite element method to consider a domain which includes an area with a complicated geometry that is difficult to approximate with quadrangles. In this work, two types of matching condition are considered which are the point-wise and the integral matching condition. The study of their convergence properties shows that the integral matching condition is the most accurate while both approaches have comparable computational costs. The matching condition can also involve a Lagrange multiplier like in (Le [START_REF] Tallec | Domain decomposition with nonmatching grids: augmented Lagrangian approach[END_REF]. For instance, Wieners and Wohlmuth [START_REF] Wieners | The coupling of mixed and conforming finite element discretizations[END_REF] proposed the coupling of Raviart-Thomas finite elements with conforming finite elements using mortar elements at the interface between subdomains with a piece-wise constant Lagrange multiplier.

Regarding the solution of wave problems with mortar elements at the interface of the subdomains, the second-order acoustic wave equation is addressed in [START_REF] Lahaye | Hybrid finite element-spectral element approximation of wave propagation problems[END_REF].

Then a variant of the classic mortar element method with Lagrange multiplier has been developed in [START_REF] Belgacem | The mortar finite element method with Lagrange multipliers[END_REF] for the Maxwell system. It is worth noting that this variant is always nonconforming because of the relaxation of matching constraints at the vertices (or the edges in 3D) of subdomains. Mortar finite elements have been also developed for solving complex elastodynamic problems in [START_REF] Casadei | A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems[END_REF]. Here SEm is coupled with a standard finite element method. The robustness and efficiency of the mortar element method is demonstrated in different cases including complex geometries approximated with finite elements. However, the coupling with mortar finite elements is accurate if the degrees of freedom of the two numerical methods are coupled at the interface. This constraint prevents from applying a full explicit timemarching scheme which is usually preferred when dealing with large scale domain. It is worth noting that in this work the size of the interface is significantly smaller than the dimension of the global problem. Nevertheless in the case of geophysical problems, it might hamper the efficiency of this approach. The use of numerical fluxes does avoid this technical step of the implementation and this point is clearly in favor of DG approximations, which naturally involve numerical fluxes. Basically, the coupling of a discontinuous method with a continuous method should prevent from this issue. For instance, it has been shown in [START_REF] Parmantier | Application of a hybrid finite difference/finite volume method to solve an automotive EMC problem[END_REF]) that numerical fluxes can be used at the interface between subdomains to get a hybrid scheme mixing Finite Difference method (FDm) with SEm for solving electromagnetic problems accurately. This work has been extended to the three-dimensional case in [START_REF] He | Hybrid finite difference/finite volume method for 3-D conducting media problems[END_REF]. Very recently, finite differences based upon Yee scheme have been coupled [START_REF] Deymier | Etude d'une méthode d'éléments finis d'ordre élevé et son hybridation avec d'autres méthodes numériques pour la simulation électromagnétique instationnaire dans un contexte industriel[END_REF] with finite volumes for solving Maxwell equations in complex media including for instance antennas.

As far as DGm is concerned, Perugia and Schötzau [START_REF] Perugia | On the coupling of local discontinuous Galerkin and conforming finite element methods[END_REF] have built and analyzed a coupling of DGm with FEm for solving elliptic problems. For the same type of equations, the coupling of local DGm with Raviart-Thomas finite elements has been achieved in [START_REF] Cockburn | Approximation of the velocity by coupling discontinuous Galerkin and mixed finite element methods for flow problems[END_REF]. The coupling of discontinuous Galerkin methods with mixed finite elements has also been applied to the computation of flows in porous media [START_REF] Girault | Coupling discontinuous Galerkin and mixed finite element discretizations using mortar finite elements[END_REF]. Then, a unified framework has been proposed in [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF] for hybridization of finite element methods in the case of second-order elliptic equations.

In this chapter, we develop a method that couples DGm and SEm. The method is described in details in the case of the elastic wave equation and it is just presented for the acoustic case. It is worth noting that the coupling involves the mixed-primal formulation of DGm. Indeed, the use of this formulation facilitates the matching with the SE formulation that only involves numerical fluxes at the interface between the two subdomains. This chapter is ending with the case of the elasto-acoustic system for which we have to consider both the numerical and the physical coupling.

DG-SEm variational formulation in the elastic case

We begin with developing the DG-SEm formulation for solving the elastic wave equation which is set in the bounded domain Ω s with boundary Γ s . We impose the boundary condition σn = 0 where n is the unitary normal vector outwardly directed to Ω s,h and defined on Γ s . We assume that Ω s is divided into two subdomains Ω 1 and Ω 2 with no overlapping and the interface separating Ω 1 from Ω 2 is denoted by Γ 1,2 . The domain is paved with a set of elements K which are composed of Ω 1 and a set of triangles (tetrahedra in 3D) covering Ω 2,h (see Figure 4.2). Ω 1 is then considered as a macro-element of the triangulation in which DGm is applied. It is matched with the elements in the triangulation of Ω 2 through numerical fluxes which are defined at the interface Γ 1,2 only. We are thus taking advantage of DG formulation, which consists in reconstructing the numerical solution from local solutions that communicate with each other through fluxes defined at the internal boundaries of the DG mesh. The coupling is then implemented by paving the macro-element with quadrangles and applying a SEm inside. Regarding the mesh denoted by T h , it is composed of cells K that are unstructured triangle/tetraedral cells and for the time being, Ω 1,h considered as one structured quadrangle cell. We keep going by writing the global DG variational formulation. Let v h ∈ V h and σ h ∈ Σ h be the global unknowns in Ω s,h . The discontinuous finite element spaces V h and Σ h are defined as:

V h = {v h ∈ L 2 (Ω s,h ) d | (v h ) |K ∈ V K (K), ∀K ∈ T h }, Σ h = {σ h ∈ L 2 (Ω s,h ) d 2 | (σ h ) |K ∈ Σ K (K), ∀K ∈ T h and (σ h ) ij = (σ h ) ji },
where V K and Σ K are two finite subspaces of (H 1 (K)) d and (H 1 (K)) d 2 .

We consider a test-function w h ∈ V h and a test-tensor ξ h ∈ Σ h . The mixed-primal DG variational formulation reads:

             Ω s,h ρ∂ t v h • w h = - Ω s,h σ h : ∇w h + Γ s,int {{σ h }} : [[w h ]], Ω s,h ∂ t σ h : ξ h = Ω s,h (Cξ h ) : ∇v h - Γ s,int {{Cξ h }} : [[v h ]]. (4.1) Γ 1,2 Ω 2,h Ω 1,h
             Γ s,int {{σ h }} : [[w h ]] = Γ 2,int {{σ h }} : [[w h ]] + Γ 1,2 {{σ h }} : [[w h ]], Γ s,int {{Cξ h }} : [[v h ]] = Γ 2,int {{Cξ h }} : [[v h ]] + Γ 1,2 {{Cξ h }} : [[v h ]]. (4.2)
Then, we plug the relations (4.2) in the formulation (4.1):

             Ω s,h ρ∂ t v h • w h = - Ω s,h σ h : ∇w h + Γ 2,int {{σ h }} : [[w h ]] + Γ 1,2 {{σ h }} : [[w h ]], Ω s,h ∂ t σ h : ξ h = Ω s,h (Cξ h ) : ∇v h - Γ 2,int {{Cξ h }} : [[v h ]] - Γ 1,2 {{Cξ h }} : [[v h ]]. (4.3)
The coupling is then implemented by considering that Ω s,h = Ω 1,h ∪ Ω 2,h . The terms on Γ 1,2 , written in blue, are the ones that couple the DG solution computed in Ω 2,h and the solution computed in Ω 1,h . The latter is obtained by applying SEm in Ω 1,h which is paved with structured quadrangles. It is worth noting that the choice of the mixed-primal DG formulation is particularly well adapted to the coupling with SEm.

DG-SEm variational formulation in the acoustic case

We introduce the spaces V h and P h :

V h = {v h ∈ L 2 (Ω f,h ) d | (v h ) |K ∈ V K (K), ∀K ∈ T h } P h = {p h ∈ H 1 (Ω f,h ) | (p h ) |K ∈ P K (K), ∀K ∈ T h }.
where P K is a subspace of H 1 (K).

Let v h and p h be the global solution in V h and P h , w h and q h in V h and P h be test functions. In addition, we apply the boundary condition v h •n = 0 on the exterior boundary of the computational domain Ω f,h . The first-order system of acoustic wave equation reads:

             ρ Ω f,h ∂ t v h • w h = - Ω f,h ∇p h • w h + Γ f,int [[p h ]] • {{w h }}, 1 c 2 0 ρ Ω f,h ∂ t p h q h = Ω f,h ∇q h • v h - Γ f,int [[q h ]] • {{v h }}.
(4.4)

By following the same approach as for the elastic wave equation, we get the following DG-SEm formulation for the first-order acoustic wave equation:

             ρ Ω f,h ∂ t v h • w h = - Ω f,h ∇p h • w h + Γ 2,int [[p h ]] • {{w h }} + Γ 1,2 [[p h ]] • {{w h }}, 1 c 2 0 ρ Ω f,h p h q h = Ω f,h ∇q h • v h - Γ 2,int [[q h ]] • {{v h }} - Γ 1,2 [[q h ]] • {{v h }}. (4.5)
It is worth noting that, in the case of a general elasto-acoustic domain, there is no reason to have a flat interface between the fluid and the solid domain. To avoid this additional difficulty, we will systematically couple DGm and SEm in the fluid area (in a subdomain of Ω f ) or in the solid area (in a subdomain of Ω s ).

Energy Study

In the following, we are going to show that the DG-SEm formulation preserves a continuous and a semi-discrete energy. Then we will construct the discrete functional based upon the Leap-Frog time scheme and we will establish that it is a conservative energy if a CFL condition is satisfied. The energy study is performed in details for the elastic formulation given in the previous section, but this study has also been carried out in the acoustic case and lead to the same conclusions.

Energy of the approximate solution

In this section we study the time stability of the approximate solution by the mean of an energy analysis. The functional is defined independently of the local approximation space and the time scheme.

Let v 1,h and σ 1,h be the discrete SEm solution in Ω 1,h and v 2,h and σ 2,h the discrete DG solution in Ω 2,h .

We introduce now v

h = v 1,h v 2,h and σ h = σ 1,h σ 2,h
the global approximate wavespeed and strain tensor. In the same way, we introduce w h and ξ h as respectively the global testfunction and test-tensor. We remind that C is a symmetrical positive tensor, which means that it can be inverted in the space of symmetrical tensor. On the exterior boundary of Ω s,h we impose the boundary condition σ h n = 0 where n is the unitary normal vector outwardly directed to the domain. Let E be the functional defined in V h and Σ h by:

E = 1 2 Ω h ρv h • v h + Ω h σ h : (C -1 )σ h .
It defines an energy since C is positive definite. We choose

w h = v h , ξ h = (C -1
)σ h in the variational formulation (4.3) and we get:

             Ω h ρ∂ t v h • v h = - Ω h σ h : ∇v h + Γ 2,int {{σ h }} : [[v h ]] + Γ 1,2 {{σ h }} : [[v h ]], Ω h ∂ t σ h : (C -1 )σ h = Ω h σ h : ∇v h - Γ 2,int {{σ h }} : [[v h ]] - Γ 1,2 {{σ h }} : [[v h ]]. (4.6)
We easily see that we have

d dt E = 0. (4.7)
The derivative in time of the energy being zero, we get that the energy of the DG-SEm scheme is constant and thus we prove that the approximate solution is stable in time.

Semi-discrete energy

We introduce the following finite dimensional spaces:

V 1,h = {v 1,h ∈ H 1 (Ω 1,h ) d | (v 1,h ) |K ∈ Q m (K) d , ∀K ∈ T 1,h } V 2,h = {v 2,h ∈ L 2 (Ω 2,h ) d | (v 2,h ) |K ∈ P m (K) d , ∀K ∈ T 2,h }, Σ 1,h = {σ 1,h ∈ L 2 (Ω 1,h ) d 2 | (σ) |K ∈ Q m (K) d 2 ∀K ∈ T 1,h , and σ ij = σ ji } Σ 2,h = {σ 2,h ∈ L 2 (Ω 2,h ) d 2 | (σ) |K ∈ P m (K) d 2 ∀K ∈ T 2,h , and σ ij = σ ji },
where Q m (K) and P m (K) are the spaces of Lagrangian polynomials of degree m defined respectively on a quadrangle (or cuboid) and on a triangle (tetrahedron).

Using the same notations as in chapter 3, the matrix system reads:

   M v h ∂ t v h + R v h σ h + R 1,2 σ h σ h = 0, M σ h ,C -1 ∂ t σ h -R σ h ,C -1 v h -R 1,2 v h ,C -1 v h = 0, (4.8a) (4.8b)
where

M v h = M v 1,h 0 0 M v 2,h , R σ h ,C -1 = R σ 1,h ,C -1 0 0 R σ 2,h ,C -1 , M σ h ,C -1 = M σ 1,h ,C -1 0 0 M σ 2,h ,C -1 , R v h = R v 1,h 0 0 R v 2,h , R 1,2 v h ,C -1 =   0 R 2,1 v 2,h ,C -1 R 1,2 v 1,h ,C -1 0   , R 1,2 σ h = 0 R 2,1 σ 1,h R 1,2 σ 2,h 0 .
The matrices M v 2,h , R v 2,h are the matrices that have been previously defined when introducing the DG discretization in Chapter 3 on page 42. Matrix M v 1,h is defined in the same way as the matrix M p h defined for SEm in Chapter 3 on page 57 and matrix R v 1,h is defined like the stiffness matrix defined for SEm on page 57. Matrices

M σ 1,h ,C -1 , M σ 2,h ,C -1 , R σ 1,h ,C -1 and R σ 2,h ,C -1
, are constructed in the same way as matrices M σ 1,h (which is built as matrix M v f,h displayed on page 57 dealing with SEm) M σ 2,h , R σ 1,h and R σ 2,h but the physical parameters are included into the mass matrices instead of the stiffness matrices.

In particular we have:

R v 1,h ,C -1 = (R σ 1,h ) T , R v 2,h ,C -1 = (R σ 2,h ) T . (4.9)
The mass matrices M σ 1,h ,C -1 and M σ 2,h ,C -1 are such that each sub-matrix M K σ 1,h ,C -1 and M K σ 2,h ,C -1 has its blocks multiplied by the corresponding coefficient of matrix C -1 :

M K σ * ,h ,C -1 = C -1 ij M K , ∀i, j = 1..d(d + 1)/2, (4.10)
where M K is defined in equation (3.24) in the section devoted to DG discretization (chapter 3 on page 42). The matrices R 1,2 v h ,C -1 and R 1,2 σ h have been introduced for the coupling and we have

(R 1,2 σ 1,h ) ijkpq = Γ 1,2 Φ q e i ⊗ nψ K p e j ⊗ e k ∀p = 1..N m , q = 1..L and R 1,2 σ 1,h = (R 2,1 v 2,h ,C -1 ) T (R 1,2 σ 2,h ) ijkpq = Γ 1,2 ψ K q e i ⊗ nϕ K p e j ⊗ e k ∀p = 1..(m + 1) d , q = 1..N m and R 1,2 σ 2,h = (R 2,1 v 1,h ,C -1 ) T
where ϕ p is a local basis function for SEm, Φ q is a global basis function for SEm and ψ ⋆ is a DG basis function. Here we have i = x, y, z, j = x, y, z, k = x, y, z, and N m is the number of basis functions in an element defined in Chapter 3 in the DG discretization. Here, n is the unitary normal vector outwardly directed to Ω 1 and defined on boundary Γ 1,2 . In addition, we have:

R 2,1 σ h = (R 2,1 v h ,C -1 ) T . (4.12)
Let us introduce the scalar product < ., . > which is defined for a pair of vectors x and

y of R L , L ∈ N by < x, y >= L i=1 x i y i . In particular, if x = Az, A is a L × L matrix: < x, y >=< Az, y >= i ( j A ij z j )y i .
According to (4.8a) and (4.8b), we have:

< ∂ t M v h v h , v h > + < R v h σ h , v h > +< R 1,2 σ h σ h , v h > = 0, < ∂ t M σ h ,C -1 σ h , σ h > -< R σ h ,C -1 v h , σ h > -< R 1,2 v h ,C -1 v h , σ h > = 0. (4.13a) (4.13b)
Let us introduce the functional E sd defined by

E sd = 1 2 < M v h v h , v h > + < M σ h ,C -1 σ h , σ h > , (4.14)
Since the mass matrices are positive definite, E sd is a quadratic form and defines thus a semi-discrete energy. Then, by adding equation (4.13a) and (4.13b), and using relation (4.12) and (4.9), it comes (4.15) which shows that the semi-discrete energy is a constant in time. Hence the semi-discrete system is stable in time.

d dt E sd = 0,

Discrete energy and CFL number

We continue our analysis by investigating the existence of a full discrete energy, which depends on the time scheme that is applied for time integration. Here, we restrict our study to the case of Leap-Frog scheme. The discrete system is given by:

M v h v n+1 h -v n h ∆t + R v h σ n+ 1 2 h + R 1,2 σ h σ n+ 1 2 h = 0, M σ h ,C -1 σ n+ 3 2 h -σ n+ 1 2 h ∆t -R σ h ,C -1 v n+1 h -R 1,2 v h ,C -1 v n+1 h = 0. (4.16a) (4.16b)
Equation (4.16b) can be rewritten at a previous time-step:

M σ h ,C -1 σ n+ 1 2 h -σ n-1 2 h ∆t -R σ h ,C -1 v n h -R 1,2 v h ,C -1 v n h = 0. (4.17)
We introduce the functional E n d that is defined by:

E n d = 1 2 < M v h v n h , v n h > + 1 2 < M σ h ,C -1 σ n-1 2 h , σ n+ 1 2 h > . (4.18)
For the time being, we can not claim that E n d is an energy because there is not absolute certainty that E n d is positive. This is due to the second term in (4.18) which involves σ h at two different time-steps. We then follow a two-stage approach: (1) evaluate

E n+1 d -E n d ; (2) prove that E n d is an energy. Step (1) : evaluation of E n+1 d -E n d . We have: E n+1 d -E n d = 1 2 < M v h v n+1 h , v n+1 h > + 1 2 < M σ h ,C -1 σ n+ 1 2 h , σ n+ 3 2 h >- 1 2 < M v h v n h , v n h > - 1 2 < M σ h ,C -1 σ n-1 2 h , σ n+ 1 2 h > Then, < M v h v n+1 h , v n+1 h > -< M v h v n h , v n h >=< M v h v n+1 h , v n+1 h -v n h > + < M v h v n+1 h , v n h > -< M v h v n h , v n h -v n+1 h > -< M v h v n h , v n+1 h > which simplifies to: < M v h v n+1 h , v n+1 h > -< M v h v n h , v n h >=< M v h (v n+1 h + v n h ), v n+1 h -v n h > or equivalently to < M v h v n+1 h , v n+1 h > -< M v h v n h , v n h >=< M v h (v n+1 h -v n h ), v n+1 h + v n h >, (4.19)
since M v h is symmetrical. We also have:

< M σ h ,C -1 σ n+ 1 2 h , σ n+ 3 2 h > -< M σ h ,C -1 σ n-1 2 h , σ n+ 1 2 h >=< M σ h ,C -1 σ n+ 1 2 h , σ n+ 3 2 h -σ n-1 2 h >, since M σ h ,C -1 is symmetrical.
Then, according to the relation :

σ n+ 3 2 h -σ n-1 2 h = (σ n+ 3 2 h -σ n+ 1 2 h ) + (σ n+ 1 2 h -σ n-1 2 h ) we get < M σ h ,C -1 σ n+ 1 2 h , σ n+ 3 2 h > -< M σ h ,C -1 σ n-1 2 h , σ n+ 1 2 h >= < M σ h ,C -1 σ n+ 1 2 h , σ n+ 3 2 h -σ n+ 1 2 h > + < M σ h ,C -1 σ n+ 1 2 h , σ n+ 1 2 h -σ n-1 2 h >
By using again the fact that M σ h ,C -1 is symmetrical, we obtain :

< M σ h ,C -1 σ n+ 1 2 h , σ n+ 3 2 h > -< M σ h ,C -1 σ n-1 2 h , σ n+ 1 2 h >= < M σ h ,C -1 (σ n+ 3 2 h -σ n+ 1 2 h ), σ n+ 1 2 h > + < M σ h ,C -1 (σ n+ 1 2 h -σ n-1 2 h ), σ n+ 1 2 h > (4.20)
We now plug the equation (4.16a) in (4.19) and (4.16b) in (4.20). First, we get

< M v h v n+1 h , v n+1 h > -< M v h v n h , v n h >=< -∆t(R v h + R 1,2 σ h )σ n+ 1 2 h , v n+1 h + v n h >, (4.21)
and we also have

< M σ h ,C -1 σ n+ 1 2 h , σ n+ 3 2 h > -< M σ h ,C -1 σ n-1 2 h , σ n+ 1 2 h >= < ∆t(R σ h ,C -1 + R 1,2 v h ,C -1 )v n+1 h , σ n+ 1 2 h > + < ∆t(R σ h ,C -1 + R 1,2 v h ,C -1 )v n h , σ n+ 1 2 h > . (4.22)
To conclude, we use (4.9) and (4.12) and we sum (4.21) and (4.22) which implies that:

E n+1 d -E n d = 0 (4.23)
Hence, the functional E n d does not vary at each iteration.

Step (2): Study the sign of E n d . For that purpose, we develop the second term in the definition of E n d , following the principle of quadrature expansions:

< M σ h ,C -1 σ n+ 1 2 h , σ n-1 2 h >=< M σ h ,C -1 σ n+ 1 2 h + σ n-1 2 h 2 , σ n+ 1 2 h + σ n-1 2 h 2 > - < M σ h ,C -1 σ n+ 1 2 h -σ n-1 2 h 2 , σ n+ 1 2 h -σ n-1 2 h 2 >
Then the term under study decomposes itself into a term which is positive since M σ h ,C -1 is positive definite and a term that we can modify thanks to (4.17). We then have

< M σ h ,C -1 σ n+ 1 2 h -σ n-1 2 h 2 , σ n+ 1 2 h -σ n-1 2 h 2 >= < ∆t 2 (R σ h ,C -1 + R 1,2 v h ,C -1 )v h , ∆t 2 M -1 σ h ,C -1 (R σ h ,C -1 + R 2,1 v h ,C -1 )v h > .
We then obtain for E n d :

E n d = 1 2 < M v h v n h , v n h > + ∆t 2 4 < (R σ h ,C -1 + R 1,2 v h ,C -1 )v n h , M -1 σ h ,C -1 (R σ h ,C -1 + R 2,1 v h ,C -1 )v h > + < M σ h ,C -1 σ n+ 1 2 h + σ n-1 2 h 2 , σ n+ 1 2 h + σ n-1 2 h 2 >.
From this last expression, we see that

E n d is positive if < M v h v n h , v n h > - ∆t 2 2 < (R σ h ,C -1 + R 1,2 v h ,C -1 )v n h , M -1 σ h ,C -1 (R σ h ,C -1 + R 2,1 v h ,C -1 )v h > (4.24)
is positive. Given that M v h is positive definite, (4.24) can be rewritten as:

< v n h , I - ∆t 2 2 M -1 2 v h (R σ h ,C -1 + R 1,2 v h ) T M -1 σ h ,C -1 (R σ h ,C -1 + R 2,1 v h ,C -1 )M -1 2 v h v n h >(4.25)
where I is the identity matrix. Let B be the matrix defined by

B = ∆t 2 2 M -1 2 v h (R σ h ,C -1 + R 1,2 v h ,C -1 ) T M -1 σ h ,C -1 (R σ h ,C -1 + R 2,1 v h ,C -1 )M -1 2 v h
Then (4.24) is positive if the matrix A = I -B is positive. We can say that B is a symmetrical matrix; indeed,

M v h is symmetrical, M σ,C -1 is symmetrical as well and since (R 1,2 v h ,C -1 ) T = R 2,1 v h ,C -1 , we have that (R σ h ,C -1 + R 2,1 v h ,C -1 ) T = R T σ h ,C -1 + R 1,2 v h ,C -1 .
We can thus define the eigenvalues of B and if we denote by λ max its larger eigenvalue, we obtain that A is positive if 1 -∆t 2 2 λ max ≥ 0. We then have: Theorem 1. Let M be the matrix defined by

M = ∆t 2 2 M -1 2 v h (R σ,C -1 + R 1,2 v h ,C -1 ) T M -1 σ h ,C -1 (R σ,C -1 + R 2,1 v h ,C -1 )M -1 2 v h .
Let λ max be the largest eigenvalue. Then, if we choose ∆t such that ∆t ≤ 2 λ max , E n d is an energy associated to the DG-SEm scheme. Moreover, it is conservative.

As a conclusion, we have been able to prove that the DG-SEm scheme is stable in time by considering functionals of energy. Regarding the discrete scheme obtained when using the Leap-frog scheme, we have obtained a result of stability as well providing the time step satisfies a CFL condition. We have not been able to develop the same study when applying the RK scheme. This is for a technical reason, the construction of a discrete energy being not that obvious in this case. In the last section of this chapter, we present the particular case of the elasto-acoustic equation to which the DG-SEm scheme is applied. It is worth noting that we will assume that the physical interface between the fluid and the solid coincides with the numerical interface limiting the region where SEm is applied from the region where DGm is applied.

Elasto-acoustic coupling

Regarding geophysical exploration, it is quite usual to consider domains that are composed of a layer of fluid and a layer of the subsurface. For addressing this kind of situation, we could solve the elastodynamic system and set the velocity v s to 0. This approach may be unstable [START_REF] Bossy | Evaluation ultrasonore de l'os cortical par transmission axiale: modélisation et expérimentation in vitro et in vivo[END_REF] and more importantly, its solution requires computing vectorial unknowns even in the fluid while the acoustic wave equation governs a scalar (pressure) and a vector (velocity). Hence, for stability issues and computational cost limitation, it turns out to be interesting to consider the elasto-acoustic equation. As previously explained, we use DGm in the area of the domain which requires flexibility of the mesh (basically in the subsurface) while we use SEm in the area that can be easily represented with quadrangles (basically the fluid). Here, we address the particular case where the interface between DG and SEm coincides with the interface between the fluid and the solid. We are thus in the spirit of [START_REF] Phillips | A coupling of mixed and discontinuous Galerkin finite-element methods for poroelasticity[END_REF] which addresses the numerical solution of the coupled fluid and mechanics in Biot's consolidation model of poroelasticity by using a mixed-finite element method for the flow variables and a DGm for the displacement. However, this is not a limitation to consider general cases. Indeed, if the fluid-solid interface is arbitrarily shaped, we can apply the DG-SEm in the fluid and a full DG in the remainder of the domain.

Let Ω be the propagation domain composed of a fluid Ω f and a solid Ω s . The two domains Ω f and Ω s are respectively paved with structured quadrangle (cuboid) cells and unstructured triangles (tetrahedra). The two discretized domains are denoted by Ω f,h and Ω s,h . We denote by T f,h the set of structured quadrilaterals in Ω f,h and we denote by T s,h the set of unstructured cells in Ω s,h We denote by Γ f,h,out the external boundary of the fluid domain and Γ s,h,out the external boundary of the solid domain and we define Γ f,s as the interface between the fluid and the solid. In addition we define Γ s,h,int the set of internal boundaries of Ω s,h . The elasto-acoustic system is given by:

     ρ∂ t v f = -∇p, 1 ρc 2 0 ∂ t p = -∇ • v f , (4.26a) (4.26b) ρ∂ t v s = -∇ • σ, ∂ t σ = -C(ǫ(v s )). (4.27a) (4.27b)
This system of equations is completed with initial conditions: σ(x, 0) = 0 in Ω s , v s (x, 0) = 0 in Ω s , v f (x, 0) = 0 in Ω f and p(x, 0) = 0 in Ω f . On the external boundary, we have the boundary conditions: v f • n f = 0 on Γ f,h,out and v s • n s = 0 on Γ s,h,out where n f and n s are the unitary normal vectors defined on the external boundaries Γ s,h,out , Γ f,h,out and outwardly directed to Ω f,h and Ω s,h . On the physical interface Γ f,s , we have the following transmission condition:

v f • n = v s • n, pn = -σn, (4.28)
where n is the unitary normal vector oriented from the fluid to the solid domain.

Variational formulation

The elasto-acoustic equation has been solved with SEm in [START_REF] Diaz | Approches analytiques et numériques de problèmes de transmission en propagation d'ondes en régime transitoire. application au couplage fluide-structure et aux méthodes de couches parfaitement adaptées[END_REF] and it has been shown that a primal-dual formulation helps avoiding the use of Lagrange multipliers that are a priori required for handling the transmission condition. Here we consider the classic SE variational formulation for the acoustic system and a mixed-dual DG formulation in the solid. Hence, only equation (4.26b) is integrated by part in the elastic system to keep the two formulations consistent. Let q, w f , w s and ξ be test-functions respectively in V f,h , V f,h , V s,h and Σ h with:

V f,h = {q ∈ H 1 (Ω f,h ) d | q |K ∈ V f (K) ∀K ∈ T f,h } V f,h = {v f ∈ L 2 (Ω f,h ) d | (v f ) |K ∈ V f (K) d ∀K ∈ T f,h }, V s,h = {v s ∈ L 2 (Ω s,h ) d | (v s ) |K ∈ V s (K) d ∀K ∈ T s,h }, Σ h = {σ ∈ L 2 (Ω s,h ) d 2 | (σ) |K ∈ Σ(K) d 2 and σ ij = σ ji , ∀K ∈ T s,h }, where V f (K), V f (K) d , V s (K) d and Σ(K) d 2 are finite subspaces of H 1 (K), H 1 (K) d and H 1 (K) d 2
The variational equations related to the elasto-acoustic system read:

                                                       ρ Ω f,h ∂ t v f • w f = - Ω f,h ∇p • w f , 1 ρc 2 0 Ω f,h ∂ t pq = Ω f,h ∇q • v f - Γ f,h,out (qn f ) • v f - Γ f,s (qn) • v f , ρ Ω s,h ∂ t v s • w s = - Ω s,h (∇ • σ) • w s + Γ s,h,int {{w s }} • [[σ]], Ω s,h ∂ t σ : ξ = Ω s,h (∇ • (Cξ)) • v s - Γ s,h,int [[Cξ]] • {{v s }}, - Γ s,h,out (Cξn s ) • v s - Γ f,s (Cξn) • v s . (4.29)
Here, the terms on the boundaries Γ s,h,out and Γ f,h,out vanish due to the boundary condition we impose. Then, we regroup the two equations involving the boundary integrals on Γ f,s :

1 ρc 2 0 Ω f,h ∂ t pq + Ω s,h ∂ t σ : ξ = Ω f,h ∇q • v f - Γ f,s (qn) • v f + Ω s,h (∇ • (Cξ)) • v s - Γ s,int [[Cξ]] • {{v s }} - Γ f,s (Cξn) • v s
Let us consider the terms on Γ f,s . First, we have :

(qn) • v f = q(v f • n) = q(v s • n) (4.30)
according to the first condition of (4.28). The other term on boundary Γ f,s can be rewritten as:

(Cξn) • v s = ((Cξn) • n)(v s • n) + ((Cξn) • τ )(v s • τ ), (4.31)
where τ is the tangential vector associated with the unitary normal vector n. Then we sum up the expressions obtained in (4.30) and (4.31), and we get

q(v s • n) + ((Cξn) • n)(v s • n) + ((Cξn) • τ )(v s • τ ) = 1 2 ((Cξn) • n + q)(v s • n + v f • n) + 1 2 ((Cξn) • n + q)(v s • n -v f • n)+ ((Cξn) • τ )(v s • τ ). (4.32)
This expression is simplified thanks to the first condition of (4.28) and we finally get :

q(v s • n) + ((Cξn) • n)(v s • n) + ((Cξn) • τ )(v s • τ ) = 1 2 ((Cξn) • n + q)(v s • n + v f • n) + ((Cξn) • τ )(v s • τ ) (4.33)
We plug this expression into the variational formulation (4.29) to get

                                                       ρ Ω f,h ∂ t v f • w f = - Ω f,h ∇p • w f , 1 ρc 2 0 Ω f,h ∂ t pq = Ω f,h ∇q • v f - 1 2 Γ f,s q(v f • n + v s • n), ρ Ω s,h ∂ t v s • w s = - Ω s,h (∇ • σ) • w s + Ω s,h,int {{w s }} • [[σ]], Ω s,h ∂ t σ : ξ = Ω s,h (∇ • (Cξ)) • v s - Γ s,h,int [[Cξ]] • {{v s }}, - 1 2 Γ f,s ((Cξn) • n)(v s • n + v f • n) - Γ f,s ((Cξn) • τ )(v s • τ ).
At this stage of the development, there are still two variational equations that do not involve terms on the interface Γ f,s . Hence, we propose to modify them by adding terms which vanish as long as we have the exact solution. These additional terms will play a role when moving on the approximate solution, in particular regarding the behavior of the energy associated with the problem. They are based upon the transmission conditions (4.28) which are imposed on Γ f,s . These terms are written in red in the following variational formulation that we use for implementing the DG-SEm for elasto-acoustic equations:

                                                                     ρ Ω f,h ∂ t v f • w f = - Ω f,h ∇p • w f + 1 2 Γ f,s ((σn) • n + p)(w f • n), 1 ρc 2 0 Ω f,h ∂ t pq = Ω f,h ∇q • v f - 1 2 Γ f,s q(v f • n + v s • n), ρ Ω s,h ∂ t v s • w s = - Ω s,h (∇ • σ) • w s + Ω s,h,int {{w s }} • [[σ]] + Γ f,s ((σn) • τ )(w s • τ ) + 1 2 Γ f,s ((σn) • n + p)(w s • n), Ω s,h ∂ t σ : ξ = Ω s,h (∇ • (Cξ)) • v s - Γ s,h,int [[Cξ]] • {{v s }} - 1 2 Γ f,s ((Cξn) • n)(v s • n + v f • n) - Γ f,s ((Cξn) • τ )(v s • τ ). (4.34)
It is worth noting that this formulation is restricted to the case where the elasto-acoustic interface is flat and coincides with the boundary dividing the SE area from the DG area.

To address general cases where the shape of the interface can vary, we implement the DG-SEm either in the fluid or in the solid.

Semi-discrete system

We introduce the discrete spaces:

V f,h = {v f,h ∈ L 2 (Ω f,h ) d | (v f,h ) |K ∈ Q m (K) d , ∀K ∈ T f,h }, V s,h = {v s,h ∈ L 2 (Ω s,h ) d | (v s,h ) |K ∈ P m (K) d , ∀K ∈ T s,h }, P f,h = {p h ∈ H 1 (Ω f,h ) | (p) |K ∈ Q m (K), ∀K ∈ T f,h }, Σ h = {σ h ∈ L 2 (Ω s,h ) d 2 | (σ) |K ∈ P m (K) d 2 , ∀K ∈ T s,h , and σ ij = σ ji },
where Q m (K) and P m (K) are the spaces of Lagrange polynomials of degree m defined respectively on quadrangles(cuboids) and on triangles (tetrahedra). Let v f,h , v s,h , σ h and p h be the approximate solutions of (4.34), sought in V f,h , V s,h , Σ h and P f,h . We consider the test functions w f,h , w s,h , ξ h and q h in V f,h , V s,h , Σ h and P f,h . The variational formulation reads as

                                                                     ρ Ω f,h ∂ t v f,h • w f,h = - Ω f,h ∇p h • w f,h + 1 2 Γ f,s ((σ h n) • n + p h )(w f,h • n), 1 ρc 2 0 Ω f,h ∂ t p h q h = Ω f,h ∇q h • v f,h - 1 2 Γ f,s q h (v f,h • n + v s,h • n), ρ Ω s,h ∂ t v s,h • w s,h = - Ω s,h (∇ • σ h ) • w s,h + Γ s,int {{w s,h }} • [[σ h ]] + Γ f,s ((σ h n) • τ )(v s,h • τ ) + 1 2 Γ f,s ((σ h n) • n + p h )(w s,h • n), Ω s,h ∂ t σ h : ξ h = Ω s,h (∇ • (Cξ h )) • v s,h - Γ s,int [[Cξ h ]] • {{v s,h }} - 1 2 Γ f,s ((Cξ h n) • n)(v s,h • n + v f,h • n) - Γ f,s ((Cξ h n) • τ )(v s,h • τ ). (4.35)
Then, by taking over the same basis functions (Φ q , ϕ p , ψ l ) as in chapter 3, we end up with the elasto-acoustic semi-discrete system:

                 ρM v f,h ∂ t V f,h + R v f,h P f,h + R f,s,1 v f,h σ h + R f,s,2 v f,h P f,h = 0, M p h ∂ t P f,h -R p h V f,h -R f,s,1 p h V f,h -R f,s,2 p h v s,h = 0, M v s,h ∂ t v s,h + R v s,h σ h + R f,s,1 v s,h σ h + R f,s,2 v s,h P f,h = 0, M σ h ∂ t σ h -R σ h v s,h -R f,s,1 σ h v s,h -R f,s,2 σ h V f,h = 0.
(4.36)

The matrices

M v f,h , R v f,h , M p h , R p h , M v s,h , R v s,h , M σ h , R σ h have been introduced in Chapter 3. The matrices R f,s,1 ⋆ and R f,s,2 ⋆
are the matrices handling both the elastoacoustic and the DG-SEm coupling. R f,s,2 v f,h and R f,s,1 p h are transposed from each other and we have:

R f,s,2 v f,h ,ipq = 1 2 Γ f,s Φ q ϕ p e i • n, ∀p = 1, .., (m + 1) d , ∀q = 1, .., L, i = x, y, z,
where L is the number of basis functions Φ

i . R f,s,2 v s,h is the transposed matrix of R f,s,2 p h , and the entries of R f,s,2 p h read as R f,s,2 p h ,ipq = 1 2 Γ f,s Φ q ψ p e i • n, ∀p = 1, .., N m , ∀q = 1, .., L i = x, y, z,
where N m is the number of basis function ψ i . The matrices R f,s,1 v f,h and R f,s,1 v s,h are represented by the following entries:

R f,s,1 vf,h,ijkpq = 1 2 Γ f,s (ψ q (e i ⊗ e j )n) • n(ϕ p e k • n) ∀p = 1, .., (m + 1) d , ∀q = 1, .., N m R f,s,1 v s,h ,ijkpq = 1 2 Γ f,s (ψ q (e i ⊗ e j )n) • n(ψ p e k • n) + (ψ q (e i ⊗ e j )n) • τ (ψ p e k • τ ), ∀p, q = 1, .., N m , with i = x, y, z, j = x, y, z, k = x, y, z. The last matrices R f,s,1 σ h , R f,s,2
σ h are dense matrices and such that:

R f,s,1 σ h = C(R f,s,1 v f,h ) T and R f,s,2 σ h = C(R f,s,1 v s,h ) T ,
where C is the matrix representation of the elasticity tensor C defined in Chapter 2.

Energy study

In order to define the energy functional associated with the system, we deal with the DG-SEm formulation in which we use the approximate solution as test function. We thus choose

w f,h = v f,h , w s,h = v s,h , ξ h = (C -1
)σ h and q h = p h in variational formulation (4.35) and we get

                                                                     ρ Ω f,h ∂ t v f,h • v f,h = - Ω f,h ∇p h • v f,h + 1 2 Γ f,s ((σ h n) • n + p h )(v f,h • n), 1 ρc 2 0 Ω f,h ∂ t p h p h = Ω f,h ∇p h • v f,h - 1 2 Γ f,s p h (v f,h • n + v s,h • n), ρ Ω s,h ∂ t v s,h • v s,h = - Ω s,h (∇ • σ h ) • v s,h + Ω s,int {{v s,h }} • [[σ h ]] + Γ f,s ((σ h n) • τ )(v s,h • τ ) + 1 2 Γ f,s ((σ h n) • n + p h )(v s,h • n), Ω s,h ∂ t σ h : σ h = Ω s,h (∇ • σ h )) • v s,h - Γ s,int [[σ h ]] • {{v s,h }} - 1 2 Γ f,s ((σ h n) • n)(v s,h • n + v f,h • n) - Γ f,s ((σ h n) • τ )(v s,h • τ ).
(4.37)

Then, we introduce

E = 1 2 ρ Ω f,h v f,h • v f,h + 1 ρc 2 0 Ω f,h p h p h + ρ Ω s,h v s,h • v s,h + Ω s,h σ h : σ h . (4.38)
By adding the four equations in (4.37), we obtain:

d dt E = 0, which shows that E is constant in time.
It is worth noting that:

• If we take a closer look to E, we can see that it can be decomposed as

E = E a + E e ,
where E a denotes the energy associated with the acoustic system defined by

E a = ρ Ω f,h v f,h • v f,h + 1 ρc 2 0 Ω f,h p h p h
and E e denotes the energy of the elastic system defined by:

E e = ρ Ω s,h v s,h • v s,h + Ω s,h σ h : σ h
We thus obtain that the energy associated with the elasto-acoustic system is defined as the addition of the acoustic energy with the elastic energy.

• If we had not added the red terms in the variational formulation, the terms on the boundary Γ f,s would not have vanished and the variation of energy would have become

d dt E = - 1 2 Γ f,s p h (v f,h • n + v s,h • n) - 1 2 Γ f,s ((σ h n) • n)(v s,h • n + v f,h • n) - Γ f,s ((σ h n) • τ )(v s,h • τ ), (4.39)
which means that we can neither guarantee energy conservation nor the stability of the scheme. Indeed, we cannot guarantee the negativity of the right hand side, and the absence of explosions or instabilities. This is why the addition of the terms in red is necessary in order to make the remaining terms vanish and obtain a system which preserves the energy.

Semi-discrete energy

Now, we look at the semi-discrete energy. By taking back the notations of the matrix system (4.36), the matrix system associated to formulation (4.37) reads as:

ρM v f,h ∂ t V f,h + R v f,h P f,h -R f,s,1 v f,h σ h -R f,s,2 v f,h P f,h = 0, (4.40a) M p h ∂ t P f,h -R p h V f,h + R f,s,1 p h V f,h + R f,s,2 p h v s,h = 0, (4.40b) M v s,h ∂ t v s,h + R v s,h σ h -R f,s,1 v s,h σ h -R f,s,2 v s,h P f,h = 0, (4.40c) M σ h ,C -1 ∂ t σ h -R v s,h ,C -1 v s,h + R f,s,1 σ h ,C -1 v s,h + R f,s,2 σ h ,C -1 V f,h = 0. (4.40d)
In equation (4.40d), the physical parameters are factorized in the mass matrix, which results in the following properties:

R v s,h ,C -1 = (R v s,h σ h ) T , R f,s,1 σ h = (R f,s,1 v f,h ) T , and R f,s,2 σ h = (R f,s,1 v s,h ) T .
We introduce the semi-discrete energy defined by

E sd,e-a =< ρM v f,h V f,h , V f,h > + < M p h P f,h , P f,h > + < M v s,h v s,h , v s,f > + < M σ h ,C -1 σ h , σ h >
Then as in the previous subsection with the approximate energy, we see that:

E sd,e-a = E sd,a + E sd,e (4.41)
Hence, according to the results obtained for the elastic and the acoustic systems apart, we have that

d dt E sd,e-a = 0.
Hence we also get stability results when considering the DG-SEm solution to the elastoacoustic system. Nevertheless, it is worth noting that this properties have been established when suing the Leap-Frog scheme. In practice, we apply a Runge-Kutta scheme in time and more precisely, we employ RK4 formula of integration. Indeed, it has been observed in [START_REF] Deriaz | Stability conditions for the numerical solution of convection-dominated problems with skew-symmetric discretizations[END_REF] that RK2 integration can be unstable. The energy analysis should thus be completed with the study of the energy functional related to RK4 scheme. Here the construction of the semi-discrete energy is still an ongoing work.

Chapter 5

DG vs DG/SEm comparison on hybrid meshes

In this chapter, we provide some comparisons between numerical solutions using DG-SEm coupling and DGm in hybrid meshes composed of triangles and structured quadrilaterals in 2D and tetrahedra and cuboids in 3D. We mainly base our comparisons on the value of the computational time. The accuracy of the numerical method is also evaluated in terms of relative error that can be computed thanks to the software Gar6more which computes analytical solutions in homogeneous and bi-layered domains. This software provides exact solutions in homogeneous and bi-layered media.

We benchmark the DG-SEm coupling with a full DG discretization in domains paved with hybrid cells composed of unstructured triangles and structured quadrilaterals. First, we will address the case of homogeneous acoustic and elastic domains. Then, we will consider a domain including an underwater salt dome in order to show the efficiency of the coupling for solving wave equations in domains having highly contrasted physical parameters. Most of the simulations are performed in 2D but we also provide a result for the elasto-acoustic wave equation in 3D. All the simulations presented here have been achieved in parallel using OpenMP with eight threads. All along this chapter and as in the previous ones, we denote by v f and p the fields computed in the acoustic domain, and by v s and σ the ones in the elastic domain. In addition, n f and n s denote the unitary normal vector outwardly directed to the fluid and solid respectively .

2D homogeneous cases

The computational domain is a square with length of 3000 meters. We consider first the case of an acoustic domain and then we address the case of a domain where elastic waves propagate. The corresponding physical parameters are displayed in Table 5.1. The waves are generated by a source term which is a second order Ricker point source located at the middle of the domain (see Figure 5.1 ) and defined by:

S = (1 -2π 2 f 2 peak (t -t peak ) 2 )e -π 2 f 2 peak (t-t peak ) 2
with f peak equals to 10Hz and t peak equals to 0.12s. The computations are based upon a hybrid mesh composed of 5250 structured quadrangles on the top and 17930 unstructured triangles in the bottom area (see Figure 5.2).

SEm and DG both employ elements of order 3. We compare the results of DG-SEm with those of a full DGm of order 3 which has been implemented in the same mesh. By this way, we compare two schemes that are based upon a hybridized discretization of the domain distributed in the same way. The simulations are done using Runge-Kutta 4 time discretization with a time step ∆t equals to 1e -4 s and the simulation duration is T = 2s. The seismograms of these numerical solutions are recorded thanks to a receiver placed at point (1000,2000) and compared with an analytical solution computed with Gar6more and recorded at the same receiver. The initial condition are : v f (x, 0) = 0 and p(x, 0) = 0 for the acoustic variables and v s (x, 0) = 0, σ(x, 0) = 0 for the elastic ones. In terms of boundary conditions, we apply v f • n f = 0 in the acoustic case and σn s = 0 in the elastic case. We observe that both simulations provide the same results. In particular, it is interesting to observe that the matching between quadrangles and triangles is transparent for the numerical waves and the change of numerical method does not affect their propagation in the sense that there are no

V p (m.s -1 ) V s (m.s -1 ) ρ(kg.m -3 ) Acoustic case 2000 0 1 Elastic case 2000 1155 2
Table 5.1: Physical parameters used of 2D homogeneous domains artificial reflections at the numerical interface separating SE solution from DG solution. The seismograms in Figure 5.5 and 5.6 show that the two numerical solutions are very close to the analytical solution. Table 5.2 shows that the relative L 2 -error is good for both cases and the computational time is in favor of DG-SEm.

Relative error (%) CPU-time(s) DGm 2.488e-3 854 DG-SEm 2.517e-3 524

Table 5.2: DGm vs DG-SEm in an acoustic domain. The full elastic case is displayed in Figure 5.7 for DG-SEm and Figure 5.8 for DGm. As for the acoustic case, we do not observe any artefact at the numerical interface between triangles and quadrangles. This is particularly important in the case of DG-SEm scheme where the communication between fluxes is handled in the variational formulation thanks to additional terms involving numerical fluxes. We have an illustration of the perfect matching between both methods. The relative L 2 errors, presented in Table 5.3, shows clearly the good accuracy of the numerical solution and as in the acoustic case, the computational time is still in favor of DG-SEm. The previous simulations allow us to validate the implementation of DG-SEm in hybrid meshes, both in the 2D acoustic and 2D elastic case. As expected, the DG-SEm is still faster than the full DGm. In the examples we have provided, the computational times are reduced by a factor greater than 1.5. To complete our study, we consider the case of a heterogeneous domain in the next section.

Stratified heterogeneous elasto-acoustic domain

Now that we have considered homogeneous test cases to investigate simulations in full acoustic and elastic propagation domains, we address the case of an heterogeneous domain composed of possibly sloped layers of solid and a water layer at the top. We are solving the elasto-acoustic wave equations and besides the interface fluid-solid, the problem involves a numerical interface which defines the transition between the application of SEm and DGm. There, the numerical interface coincides with the physical interface. The domain of interest is a 4000x5000 rectangle and is depicted on Figure 5.9, the physical parameters associated to each area being displayed in Table 5.4. This domain is meshed using a hybrid mesh composed of 47514 unstructured triangles cells all located on the solid part and 1600 squares with 800 squares in the water layer and 800 squares in a layer located in the bottom of the solid domain (see Figure 5.10). We use a low number of quadrangles because we aim at testing the coupling with SEm of high-order. In this example, we will also test the use of DG-SEm in an elastic medium which corresponds to the lower part of the computational domain.

In DG-SEm coupling, DGm is employed at order 2 and SEm is applied at order 8. The reference solution is obtained with a full DGm at order 8. case, we do not dispose of a solution computed with Gar6more. We have also computed the wave fields by applying a full DGm in the hybrid mesh, that is DGm at order 2 in the triangle area and DGm at order 8 in the square area. We use a RK4 time-scheme with a time-step equals to 1e -4 s, for a simulation done during 2s. We consider the following initial conditions: v s (x, 0) = 0 and σ(x, 0) = 0 in the solid domain and v f (x, 0) = 0 and p(x, 0) = 0 in the fluid. In terms of boundary conditions we have v f • n = 0 on the external boundary of the fluid and v s • n = 0 on the external boundary of the solid. The solution is recorded at a receiver placed at point (840, 4100) and we will compare it to the reference solution obtained with the full DG discretization at order 8. The propagation phenomenon starts with a second order Ricker source defined the same way as the one used in homogeneous cases, with a central frequency f peak equals to 10Hz and a delay t peak equals 0.12s. It is placed at point (2000,4250). The numerical results are depicted in Figure 5.11 for the DG-SEm and Figure 5.12 for the DGm. Table 5.5 presents the relative L 2 -error and the computational time. We observe that in this case we do not have a significant gain of computational time and the error is in favor of DGm. Regarding the values of CPU times, they might be surprising regarding the ones obtained in homogeneous domains which are always clearly in favor of DG-SEm.

Here they are of the same order. This might be explained by the fact that we use a SEm of order 8 applied in areas paved with coarse cells. Hence the number of quadrangles is low as compared with the number of triangles. Moreover, the order of approximation is high as compared to the one of DGm. It would be interesting to relaunch this simulation by having a higher number of quadrangles while decreasing the order of approximation in the SEm. We have relaunched the simulation in the same mesh but using this time DGm at order 3 and SEm at order 5 for the coupling and DGm at order 3 in the triangles and order 5 Table 5.5: Heterogeneous domain, DGm of order 2 in triangles and 8 in squares, DG-SEm with DG at order 2 and SEm at order 8, RK4 scheme in time in the squares when using full DG discretization. Table 5.6 presents the L 2error and computational times for this case. We observe that once again we do not have significant gain of computational time using DG-SEm instead of DGm. In addition, we observe an increasing in the L 2 -error using DGm which is logical because we use a lower order than before, but the error of DG-SEm has improved. The idea that SEm is at its best when employed at order 5 has already been illustrated in [START_REF] Durufle | Intégration numérique et éléments finis d'ordre élevé appliqués aux équations de Maxwell en régime harmonique[END_REF]. Here, we see that SEm at order 5 performs as well as SEm at order 8. But DG-SEm and DGm have still the same CPU time and this should be due to the fact that the mesh contains a low number of quadrangles. It would be interesting to redo these computations by increasing the number of quadrangles.

Relative error (%) CPU-time(s) DGm 1.216e-2 2757 DG-SEm 5.131e-2 2818

Table 5.6: Heterogeneous domain, DGm of order 3 in triangles and 5 in squares, DG-SEm with DG at order 3 and SEm at order 5, RK4 scheme in time

2D domain including a salt dome

Here we consider a synthetic domain that the team uses to benchmark. It has been constructed by L.Boillot [START_REF] Boillot | Absorbing boundary conditions for 3D elastic TTI modeling[END_REF] and offers the advantage of including Tilted Transverse Isotropy (TTI) anisotropy in addition to heterogeneities. Hence, we solve the following equations:

     ρ∂ t v f = -∇p 1 ρc 2 0 ∂ t p = -∇ • v f
(5.1a)

(5.1b)

ρ∂ t v s = -∇ • σ ∂ t σ = -C(ǫ(v s ))
(5.2a)

(5.2b)

where the tensor C in the case of TTI anisotropy is given by

C T T I =          C 11 C 12 C 13 C 14 C 15 C 16 C 22 C 23 C 24 C 25 C 26 C 33 C 34 C 35 C 36 C 44 C 45 C 46 C 55 C 56 C 66          .
As explained in Chapter 2 the coefficients of this matrix can be computed from the parameter values displayed in Table 5.7. The initial conditions are v f (x, 0) = 0 in Ω f , v s (x, 0) = 0 in Ω s , p(x, 0) = 0 in Ω f , and σ(x, 0) = 0 in Ω s . We apply the boundary conditions v f • n f = 0 on Γ f × (0, T ) and v s • n s = 0 on Γ s × (0, T ).

In Figure 5.13, we provide a picture of the domain, the physical and anisotropic parameters of each area being displayed in Table 5.7. According to the characteristics of each area of the domain, we use a mesh which is composed of:

• structured quadrangle cells in the water layer, with the exception of a thin layer limited by the interface with water sand; there is no topography to take into account, the use of quadrangles is thus adapted.

• unstructured triangles in the thin water layer touching the upper boundary of the water sand; this area forms the transition zone between quadrangles and triangles. Its upper boundary is the place where SEm and DGm communicate through appropriate fluxes.

• unstructured triangles in the remainder of the domain.

The mesh is composed of 53969 triangles and 21000 squares and is depicted in Figure 5.14.

The point source is a Ricker function which is placed on the water layer at point (1500,2650). This function is equal to

S = (1 -2π 2 f 2 peak (t -t peak ) 2 )e -π 2 f 2 peak (t-t peak ) 2
where f peak is equal to 10Hz and t peak = 0.12s. We still use RK4 scheme in time with a ∆t equal to 3e -5 and the simulation duration is T = 2s. The reference solution is computed with Elasticus using a full DG scheme of order 5. The DG-SEm coupling involves a SEm and a DGm which are both of order 3. We compare the DG-SEm solution with a DGm solution obtained with a full DGm scheme of order 3. Regarding the seismograms, they are obtained from recordings performed at a receiver placed in the water layer at point (1000, 2500). The positions of the source and the receiver are depicted in Figure 5.15. 5.17 presents the result using full DGm. We see that we get the same results for both simulations. In particular, the combination of the two different meshes does not generate any numerical artefact. The domain is meshed with 50000 cuboids in the water layer and 96511 unstructured tetrahedra in the elastic part. Figure 5.21 pictures a cross-section of the mesh used for the simulation. The slice shows unstructured tetrahedra that are used in the inferior area of the domain. The slice also highlights the transition zone matching structured hexahedra with unstructured tetrahedra. For the numerical simulation, we use a second order Ricker source:

S = (1 -2π 2 f 2 peak (t -t peak ) 2 )e -π 2 f 2 peak (t-t peak ) 2
with a peak frequency equal to f peak = 5Hz and a delay in time t peak equal to 0.24s. This point source is placed in the water layer at point (500,500,1150). The initial conditions are the following:

v s (x, 0) = 0 in Ω s , v f (x, 0) = 0 in Ω f , p(x, 0) = 0 in Ω f , σ(x, 0) = 0 in Ω s . The boundary conditions are: v f • n f = 0 on Γ f × (0, T ) and v s • n s = 0 on Γ s × (0, T ).
The time discretization is carried out with RK4 scheme with a time-step ∆t = 1e -5.

The final time of simulation is T = 0.6s.

Figure 5.22 pictures the numerical results. As in the two-dimensional case, we can see that the matching between the two computational subdomains is perfect. We can see that the propagation is not polluted by artefacts appearing in particular at the fluidsolid interface. Moreover, some more numerical experiments are necessary to validate this simulation, for instance, it would be interesting to have the same simulation with full DGm. This could provide a way for validating the computation and then it would be relevant to increase the orders of approximation to assess the computational performances of the coupling in 3D.

Conclusion

In this chapter, we have presented several comparisons between the DG-SEm scheme and DG used in hybrid meshes in 2D. All these comparisons lead to the same conclusion: there is a clear reduction (by a factor always larger than 2) of the computational times with DG-SEm without any significant loss of accuracy even if DGm provides a more accurate solution. In particular, we have considered a synthetic domain which sounds like geophysical cases including underwater salt dome. In this case, the computational time of DG-SEm is twice inferior to the one when we use DGm while keeping a very good level of accuracy. We have also validated the DG-SEm in an elastic domain. This example should be studied in more details to analyze if the ratio of triangles versus quadrangles can be optimized to improve the computational performances of DG-SEm. We have also observed that it could exist an optimal order for SEm beyond which the SEm performance is not improving. In 3D, we have presented a numerical result which is very promising and deserves being completed to deliver a relevant numerical investigation of the 3D case. Another important issue is the truncation of the computational domain. In the experiments that we have performed until now, we have considered boundary conditions that generate reflections and thus may pollute the numerical solution. In the following chapter, we address the question of introducing non reflecting boundary conditions.

Chapter 6

Domain truncation conditions

State of the art

The numerical simulation of seismic waves requires introducing domain truncation conditions. The question of truncating the propagation domain of waves interests researchers for a while and it turns out that an undeniable answer is still expected. From a theoretical point of view, we dispose of transparent boundary conditions (TBC) (see [START_REF] Ionescu | Transparent boundary conditions for wave propagation on unbounded domains[END_REF] and their references therein). They are exact in the sense that they allow for the computation of the exact wave field in a bounded region limited by a boundary on which the TBC is set. Unfortunately, TBCs are generally difficult to implement and more importantly, they generate extra computational costs that make the computations unaffordable. Extra computational costs are essentially due to the fact that TBS are governed by global operators both in time and space. Hence TBCs are rather approximated to get local operators requiring lower computational burden. This approach leads to the construction of absorbing boundary conditions (ABC). Many works deal with the construction of ABCs since the pioneering works of Clayton, Engquist and Majda [START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave equations[END_REF][START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF] at the end of the 1970's who introduced ABCs for acoustic and elastic waves. Then the issue of stability has been addressed by Reynolds [START_REF] Reynolds | Boundary conditions for the numerical solution of wave propagation problems[END_REF] who proposed a stable ABC of order 1 for elastic waves. It is worth noting that the construction of ABCs for elastic waves is much more difficult than for acoustic waves because the equations are vectorial and govern P and S waves with possible conversions of P into S and vice versa. But Higdon [START_REF] Higdon | Absorbing boundary conditions for elastic waves[END_REF] succeeded in constructing high order ABCs for elastic waves that have been widely used by geophysical community and are based on the combination of low order ABCs. ABCs have proven to be very efficient and suitable for arbitrarily shaped surfaces but they are not as absorbing than expected. Actually, one could think that increasing the order of the ABC makes the condition more absorbing but it turns out that the system of governing equations coupled with high-order ABC is not surely well-posed as illustrated in [START_REF] Trefethen | Well-posedness of one-way wave equations and absorbing boundary conditions[END_REF]. Moreover, their discretization does not always lead to stable numerical schemes, in particular with DG approximations [START_REF] Duprat | Développement de conditions absorbantes pour des problèmes de propagation d'ondes en milieux complexes[END_REF]. This observation motivated later on the construction of stable ABCS for elastic waves which are limited to first order approximation and the resulting boundary value problem is well-posed [START_REF] Barucq | Absorbing boundary conditions for 2D tilted transverse isotropic elastic media[END_REF].

In the light of those problems, many authors have proposed solution methodologies like in [START_REF] Cerjan | A nonreflecting boundary condition for discrete acoustic and elastic wave equations[END_REF] where it is proposed to solve one-way-like wave equations inside a strip surrounding the computational domain. The resulting ABC is based on progressive reduction of amplitude in the strip. This method is now commonly used as a window taper function. The idea of introducing a layer in which absorption acts has been discussed again by Bérenger in [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Berenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF] who had the idea of modifying the governing equations inside the external layer only. The modification amounts to glue the external layer to the original propagation domain through constitutive parameters that match perfectly with the original ones at the interface with the absorbing layer and that represent an ideal material in the layer where waves are absorbed exponentially. Bérenger introduced his method as Perfect Matched Layers (PML) and it was the beginning of many works devoted to the design of PMLs for different systems of equations, Maxwell system being the first one. Regarding elastic waves, interesting results are in [START_REF] Peng | An optimal absorbing boundary condition for elastic wave modeling[END_REF][START_REF] Hastings | Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation[END_REF] while acoustic waves have been considered in [START_REF] Liu | The perfectly matched layer for acoustic waves in absorptive media[END_REF]. These contributions have been followed by mathematical and numerical analysis like in [START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF] and in [START_REF] Hesthaven | On the analysis and construction of perfectly matched layers for the linearized euler equations[END_REF] for the linearized Euler equation. The literature is particularly rich on the implementation of PMLs in a SEm (see [START_REF] Komatitsch | A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation[END_REF][START_REF] Fauqueux | Eléments finis mixtes spectraux et couches absorbantes parfaitement adaptées pour la propagation d'ondes élastiques en régime transitoire[END_REF][START_REF] Festa | The Newmark scheme as velocity-stress time-staggering: an efficient PML implementation for spectral element simulations of elastodynamics[END_REF][START_REF] Madec | Méthode des éléments spectraux pour la propagation d'ondes sismiques en milieu géologique fluide-solide avec pas de temps locaux et couches absorbantes parfaitement adaptées C-PML[END_REF][START_REF] Martin | A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation[END_REF] and their references) and the PMLs turn out to be very efficient to absorb both elastic waves with non grazing incidence and surface waves. Nevertheless, to address the issue of grazing incidence, finite differences are considered most of the time. In this regard, Komatitsch and Martin [START_REF] Komatitsch | An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[END_REF] have proposed a new unsplit PML for the seismic wave equation which calls for a convolutional technique. The new PML turns out to perform better than classical PMLs at grazing incidence for the same memory cost. However, it is also shown that Convolutional Perfectly Matched Layers (C-PML) is unstable in anisotropic materials, as the classical PML. Later on, the C-PMLs have been adapted to be implemented in higher order time schemes [START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML)[END_REF]. Such absorbing conditions are called Auxiliary Differential Equation PML (ADE-PML). They keep the advantage of C-PMLs at grazing incidence and are based on Runge-Kutta time-stepping scheme. More importantly, it is shown that the ADE-PML formulation is numerically stable up to 100,000 time steps for isotropic media. In this chapter, we propose the application of the DG-SEm coupled with ABC or PML. This is an important question since we only validate the DG-SEm with surface-free conditions on the boundary of the computational domain. We focus first on the ABC. After introducing the mathematical formulation including an ABC, we present numerical results in homogeneous and heterogeneous domains formerly presented in chapter 5. We observe that some reflections may pollute the simulations and then we go to the implementation of a PML. We consider a classic unsplit PML and the ADE-PML proposed in [START_REF] Martin | A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation[END_REF]. We first test the PML formulations when solving the equations with a full DGm based upon a quadrangle mesh. The experiment becomes unstable with time while it is stable when applying SEm. We get this result both in a homogeneous and heterogeneous domain. Then we implement the ADE-PML in the DG-SEm and we can see, by performing a 2D and a 3D simulation that the SEm stabilizes the PML.

Absorbing boundary conditions

Absorbing Boundary Conditions (ABC) are studied for a while, since the pioneering work of [START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF]. When an analytical solution is available, their construction is quite straightforward since it is possible to exhibit an exact transmission condition which represents the transmission of a wave from a domain to one other and eventually the reflections if any. In the case of truncating the computational domain, ideal conditions are those which force the numerical solution to get out of the domain. By this way, the boundary used for truncation is artificial and does not interfere with the numerical solution. Unfortunately, ideal conditions are generally represented by infinite series which cannot be implemented as they are. ABCs are then constructed by truncating the infinite series and their accuracy depends on the order of truncation. However, in general cases (where the artificial surface is complex or the domain is heterogeneous, etc), we do not dispose of analytical solutions and the construction of ABCs is much more technical. [START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF] have shown that it is possible to obtain ABCs from microlocal analysis of the acoustic wave equations in the vicinity of the artificial boundary and their work has been followed by many authors, some of them extending their idea to a larger scope of equations.

ABC for acoustic isotropic medium in 2D

We keep adopting the notation v f for the velocity and p for the pressure in the fluid domain, v s for the velocity and σ for the stress tensor in the elastic domain. n denotes the unitary normal vector outwardly directed to the domain. The simplest ABC that is used for truncating an acoustic isotropic domain is well-known when applied to the most common formulation of the acoustic wave equation which is given in the plane (x, z) by:

1 c 2 ∂ 2 t p - 1 ρ (∂ 2 x p + ∂ 2 z p) = 0. (6.1)
In the above notations, ρ denotes the density of the fluid and c stands for the propagation velocity. Then the ABC reads:

1 c ∂ t p + ∂ n p = 0 (6.2)
where n denotes the unitary normal vector outwardly directed to the computational domain and associated with the boundary on which the ABC is set. It is easy to rewrite the ABC for the first order formulation of the acoustic wave system by plugging its equations in the formulation of the ABC. Indeed, the first order system says that:

       ρ∂ t v f,x = -∂ x p, ρ∂ t v f,z = -∂ z p, 1 ρc 2 ∂ t p = -∂ x v f,x -∂ z v f,z , (6.3)
Then by taking (6.3) into account, (6.2) modifies to:

1 c ∂ t p = ρ∂ t v f • n (6.4)
We then obtain that 1 c pρv f • n is constant over the time which implies that:

1 c p = ρv f • n (6.5)
when the initial conditions are zero. This is the condition that we use for truncating a fluid when solving the first order formulation of the acoustic wave equation.

Implementation of the ABC into DG-SEm formulation: acoustic case

In the case of the acoustic DG-SEm, we have two equations which combine the pressure p and the wavespeed v f . The approximate pressure is continuous in the whole domain while the velocity is locally continuous and discontinuous across each element of the triangulation. To keep the local (element per element) structure of the associated discrete system, we have decided to implement the ABC in the equation governed by v f which is discontinuous both for DGm and SEm. By taking back the notations of Chapter 4, we thus consider the equation:

ρ Ω f,h ∂ t v h • w h = - Ω f,h ∇p h • w h + Γ 2,int [[p h ]] • {{w h }} + Γ 1,2 [[p h ]] • {{w h }} (6.6)
This equation does not involve any term on the external boundary of the domain since it does not result from an integration by parts. We then add the ABC by imposing the condition weakly through an additional term. The equation is then replaced by:

ρ Ω f,h ∂ t v h • w h = - Ω f,h ∇p h • w h + Γ 2,int [[p h ]] • {{w h }}+ Γ 1,2 [[p h ]] • {{w h }} - Γ f,h,out ρc(v f • n f )(w h • n f ) + Γ f,h,out p h (w h • n f ) (6.7)
where Γ f,h,out stands for the exterior boundary of Ω f,h .

Implementation of the ABC into DG-SEm formulation: elastic case

We still turn back to the notations of Chapter 5. Then, Ω s,h denotes the computational domain, which is paved with elements K. We denote by Γ s,int the set of interior edges/faces and Γ s,out the exterior boundary that is introduced to truncate Ω s,h . According to [START_REF] Cohen | Computational methods for transient analysis[END_REF], the simplest ABC for isotropic elastic media reads as σn = Bv s , (6.8)

where n stands for the unitary normal vector outwardly directed to Ω s,h and defined on Γ s,out . The matrix B is symmetric, positive definite and depends on the Lamé parameters or equivalently on the velocities V P and V S . For instance in the 2D case, we have: (6.10) where τ is the tangential vector associated with n.

(σn) • n = -ρV p v s • n (6.9) (σn) • τ = -ρV s v s • τ,
As in the acoustic case, we add the ABC in the equation governing the variable that is discontinuous both in DGm and SEm. In the elastic case, it is σ h which is governed by the equation: 6.11) In 6.8, the ABC gives the σ as a function of the velocity v s . But since B is positive definite, we can define its inverse B -1 and write the ABC in the form:

Ω s,h ∂ t σ h : ξ h = Ω s,h ((Cξ h )) : ∇v h - Γ s,int Cξ h : [[v h ]], (
v s = B -1 σn
Hence, we get the modified equation including the ABC as follows:

Ω s,h ∂ t σ h : ξ h = Ω s,h ((Cξ h )) : ∇v h - Γ s,int Cξ h : [[v h ]] + Γs,out Cξ h • (B -1 σ h n -v s ),
(6.12)

Numerical results

The computational domain is a square of side 3000 meters (see Figure 6.1). 3000m 3000m • src Figure 6.1: Domain of computation We consider a homogeneous domain which is paved with cells composed of 17930 unstructured triangles and 5250 squares. The squares are located on the top of the domain (see Figure 6.2). We use a DG-SEm scheme of order 3 in the sense that both DGm and SEm are of order 3 in space. The time integration is performed with a 4-th order Runge-Kutta (RK4) time scheme with a CFL number ∆t equals to 1e -4 and the simulation duration T is T = 2s. The simulations are done using OpenMP on eight threads. Moreover, the initial conditions are v f (x, 0) = 0 and p(x, 0) = 0 for the acoustic case and v s (x, 0) = 0 and σ(x, 0) = 0 for the elastic case. The ABC is applied on each side of the domain. The source term is a Ricker function of order 2 placed at point (1500,1500) and defined as:

S = (1 -2π 2 f 2 peak ((t -t peak ) 2 )e -π 2 f 2 peak ((t-t peak ) 2
with f peak = 10Hz and t peak = 0.12s. The propagation velocity in the acoustic domain is equal to 2000m.s -1 and the density ρ is 1kg.m -3 . Figures 6.3 and 6.4 picture the results obtained at different times of simulation. The first one illustrates the use of the ABC in DG-SEm scheme while the second one is obtained with DGm and the ABC. We can see that both simulations provide the same snapshots and at first glance, we can observe that the external boundary let the waves out. Nevertheless, if we change the scale of Figure 6.2: Hybrid mesh used for the homogeneous cases.

amplitude to flash lower values as it is done at the third line of the two figures, we observe remaining reflections at grazing incidence. This phenomena is classical when using ABCs, in particular when they are of low order. We present similar results in an isotropic elastic domain where the P-wavespeed is equal to 2000m.s -1 , the S-wavespeed is equal to 1555m.s -1 and the density ρ is equal to 1kg.m -3 . Snapshots are depicted in Figure 6.5 and 6.6 where we respectively apply DG-SEm and DGm with ABC and the same conclusions hold.

Numerical results in a heterogeneous case

We consider the propagation domain which has been described in Chapter 5 and contains an underwater salt dome. The equations that are solved are the elasto-acoustic wave equations. The domain is discretized with the same mesh depicted in Chapter 5 in Figure 5.14. As in chapter 5, the interface between the areas in which DGm and SEm are applied is in the fluid domain, near the physical interface separating the fluid and the solid region. Both DGm and SEm are of order 3 in space and we use a RK4 scheme with a CFL number equals to 3e -5 s. The final time of the simulation here is T = 2s. The simulations are done using OpenMP on eight threads. In terms of initial conditions, we have: v s (x, 0) = 0, v f (x, 0) = 0, p(x, 0) = 0, σ(x, 0) = 0, and the ABC is set on each boundary. As source term, we use the same Ricker function as in the homogeneous domain with f peak also equal to 10Hz and a delay t peak equal to 0.12s. This point source is placed at point (1500,2550) in the layer of water. Snapshots are presented in Figure 6.7. We observe that the absorption in both acoustic and elastic domain is acceptable.

We have obtained exactly the same results when using DGm in place of DG-SEm. As in the homogeneous case, we have also remarked that reflections tend to pollute the numerical solution. Hence it is very interesting to dispose of Perfectly Matched Layer facilities. PMLs have indeed demonstrated a better efficiency than ABCs of low order in many examples. This is the subject of the next section

Perfectly Matched Layers (PML)

Perfectly Matched Layers (PML) have been developed to offer an efficient alternative to ABCs which are not able to absorb waves in each situation, in particular those which impinge the boundary at incidence far from the normal one. PMLs have been introduced by Bérenger [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] for the Maxwell equations and his pioneering work has been followed by many works which propose different versions of PMLs more or less difficult to implement. More importantly, they distinguish themselves from their associated computational cost and/or the stability of the computations. We have first implemented both classical PML and C-PML formulation for DGm.

Classic PML formulation

We remind here the first-order acoustic system:

           ρ ∂ ∂ t v f = -∇p, ∂ ∂ t p = -∇ • v f . (6.13)
To present the first version of PMLs, we retain the interpretation provided by F. Collino [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]. A PML formulation is based upon rewriting the above system in a coordinate system which amounts modifying the derivative with respect to the variable along which the absorption occurs. For that purpose it is convenient to work in the frequency domain on which the acoustic system reads:

       ρiωv f = -∇p, iωp = -∇ • v f . (6.14)
where ω the pulsation and i is the imaginary number satisfying i 2 = -1. To explain the construction of the PML equations, we consider the case of absorption in the x direction which means that we will modify the coordinate x only. The new system of coordinates involves a new variable denoted by x which is defined by:

x

(x) = x - i ω x 0 d x (s)ds. (6.15)
The function d x is a damping profile which satisfies d x (x) = 0 outside the PML domain and is positive inside the layer. It is meant to absorb waves propagation inside the PML. Then we have: .16) with:

∂ x = iω iω + d x ∂ x = 1 s x ∂ x . ( 6
s x = iω iω + d x = 1 + d x iω (6.17)
In the case where the PML is applied in a direction parallel to the z-axis, we define a similar change of coordinates for the z variable: .18) and we define in a similar way ∂ z and s z : (6.19) with:

z(z) = z - i ω z 0 d z (s)ds. ( 6
∂ z = iω iω + d z ∂ z = 1 s z ∂ z .
s z = iω iω + d z = 1 + d z iω (6.20)
Then, we get a PML formulation of the harmonic system by replacing all the derivatives ∂ x and ∂ z respectively by ∂ x and ∂ z :

       iωρv f,x = -∂ xp iωρv f,z = -∂ z p, iωp = -∂ xv f,x -∂ z v f,z (6.21)
and it is worth noting that (6.21) coincides with the initial problem (6.14) outside the PML. To obtain the formulation of the PML system in the time domain, we follow a two-stage approach: first we rewrite (6.21) in terms of x and z. Then we define auxiliary unknowns (v 1 , v 2 , p 1 , p 2 ) and by applying an inverse Fourier transform to go back to the time variable, we get the classical split PML formulation:

                                       (∂ t + d x )ρv 1 f,x = -∂ x p (∂ t + d z )ρv 2 f,x = -∂ z p (∂ t + d x )ρv 1 f,z = -∂ x p (∂ t + d z )ρv 2 f,z = -∂ z p (∂ t + d x )p 1 = ∂ x (v f,x + v f,z ) (∂ t + d z )p 2 = ∂ z (v f,x + v f,z ) v f,x = v 1 f,x + v 2 f,x v f,z = v 1 f,z + v 2 f,z p = p 1 + p 2 (6.22)
It has been shown in [START_REF] Komatitsch | An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[END_REF] that this classical PML formulation does not give satisfactory results when it comes to waves at grazing incidence. This can be overcome by the use of Convolutional-PML (C-PML) formulation which is presented in the next section.

Convolutional-PML (C-PML)

The C-PML formulation has been introduced by Roden and Gedney in [START_REF] Roden | Convolution PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media[END_REF] for Maxwell equations and then developed by Komatitsch and Martin in [START_REF] Komatitsch | An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[END_REF] for elastodynamic problems in strong formulation and in [START_REF] Martin | A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation[END_REF] for a weak formulation. The C-PML technique consists in modifying the choice of the stretching function s x define on equation (6.17) by introducing α x ≥ 0 and κ x ≥ 1 two real variables such as:

s x = κ x + d x α x + iω (6.23)
The above change in the stretching variable amounts in some sense to replace the pulsation by a complex number. By this way, the absorption of waves is increased by considering evanescent waves as well. As shown in [START_REF] Komatitsch | An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[END_REF], the value of the parameter κ x can be taken equal to one without any impact on the accuracy of the simulation. Hence we consider

s x = 1 + d x α x + iω (6.24)
Then if sx (t) is the inverse Fourier transform of 1 sx , we can express the quantity ∂ x as:

∂ x = sx (t) ⋆ ∂ x (6.25)
where ⋆ stands for the convolution product. After a rewriting of (6.24) as:

1 s x = 1 - d x α x + iω + 1
we can compute sx (t) and we get: sx (t) = δ(t) + d x e -(dx+αx)t H(t).

(6.26)

where δ(t) denotes the Dirac measure and H(t) is the distribution of Heaviside. Hence, we deduce that: (6.27) where ζ x (t) is given by: ζ x (t) = d x e -(dx+αx)t (6.28)

∂ x = ∂ x + ζ x (t) ⋆ ∂ x .
In practice, ζ x (t)⋆∂ x is approximated at iterate n∆t. In [START_REF] Martin | A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation[END_REF], the recursive convolution method of Luebbers and Hunsberger (see [START_REF] Luebbers | FDTD for Nth-order dispersive media[END_REF] for details) is applied. It is restricted to the case of zero initial data and reads: Besides including more absorption, one interest of C-PML is that they do not use split unknowns. However they are not straightforwardly adapted with time-stepping schemes of order higher than 2.

Auxiliary Differential Equation-PML (ADE-PML)

The idea of ADE-PML has been proposed in [START_REF] Gedney | An auxiliary differential equation formulation for the complex-frequency shifted PML[END_REF] and applied to seismic wave propagation in [START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML)[END_REF]. The idea is to remove the computation of the convolution product by reformulating the memory variable defined in the previous section in a similar form which can be interpreted as a time evolution of the memory variable.

Numerical experiments

In this section, we perform some numerical experiments to assess the performance of DG-SEm when solving the PML formulation of acoustic and elastic wave equations. It is worth noting that the layer introduced for absorbing the waves is systematically paved with quadrangles or cuboids. In practice, the layer does not contain any complex shape that could justify the use of triangles or tetrahedra. Since we aim at solving the PML formulation with DG-SEm applied on quadrangle mesh, we proceed step by step and first test the PML formulation with full DGm.

PML formulation solved with full DGm

We consider a 3000mx3000m square domain. It is a homogeneous elastic domain with V p = 2500m.s -1 , V s = 1250m.s -1 and ρ = 1kg.m -3 . This domain is discretized with 20000 quadrangles and we apply DGm of order 3. The time discretization is done with RK4 scheme with a time-step equals to 2.4e -4 and the simulation stands during T = 10 seconds. The PML is ended with the boundary condition σn = 0. The PML layer is placed at the top boundary and is 400 m thick (see Figure 6.8). On the other side of the domain, we apply a free-surface condition σn = 0. We have performed simulations with the two PML formulations: classic PML with ADE-PML and CPML with ADE-PML. As

• src 400m Figure 6.8: Position of the point source and the PML layer (in blue), homogeneous elastic domain in [START_REF] Komatitsch | An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[END_REF] and [START_REF] Madec | Méthode des éléments spectraux pour la propagation d'ondes sismiques en milieu géologique fluide-solide avec pas de temps locaux et couches absorbantes parfaitement adaptées C-PML[END_REF], the damping profile is polynomial. In the PML, the waves are absorbed in the direction of the z axis so we have d z defined by:

d z = Cd 0 (z/L) N (6.34)
where L is the thickness of the PML (here, L = 400), N is the degree of absorption. We also have d 0 = -(N + 1)V max P log(R c )/(2L) as in [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF] where R c is the theoretical coefficient of reflection, V max P is the maximal P -wavespeed in the propagation domain and C is a constant. In the case where we consider a C-PML, the polynomial α z has the same form as the one used in [START_REF] Martin | A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation[END_REF]:

α z = α max [1 -(z/L) m ] (6.35)
with α max = πf peak (see [START_REF] Komatitsch | An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[END_REF]) and f peak is the dominant frequency of the source. Here we have: L = 400m, N = 2, C = 1, V max p = 2000m.s -1 , m = 1, and R c = 0.001%

We set the following initial conditions: v s (x, 0) = 0 and σ(x, 0) = 0. The source we apply is a second order Ricker point source defined as: S = (1 -2π 2 f 2 peak (tt peak ) 2 )e -π 2 f 2 peak (t-t peak ) 2 and we take f peak = 10Hz and t peak = 0.12s. It is placed at the middle of the domain at point (1500,1500). Figure 6.9 depicts the numerical results obtained with the classical PML formulation and Figure 6.10 the ones corresponding to C-PML formulation. We observe that for both cases, instabilities appear on the PML layers which lead to an explosion of our solution. Then we move on with the same computations but this time, we use the DG-SEm method, the SEm being applied inside the PML. We use the same mesh as in Chapter 5. Both DGm and SEm are applied at order three, the discretization in time is done using RK4 scheme with ∆t = 1e -4 and the simulation lasts until T = 10s. The PML layer and the position of the source are the same as depicted in Figure 6.8. Moreover, to be consistent we take the same boundary conditions, the same source and the same PML and C-PML parameters. In Figure 6.11 we present the numerical results using classical PML formalism. We can see that some instabilities appear after 7s of simulation. On the other hand, we can see in Figure 6.12 than when we use the C-PML formalism, the wave is correctly absorbed and there are no instabilities appearing in the PML layer during the whole simulation. We have the confirmation as formerly observed in [START_REF] Martin | A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation[END_REF] that the addition of α x in the PML formulation makes the PML stable and more efficient at grazing incidence. Now the difference between the two series of simulations is the use of SEm inside the PML. This numerical experiment is thus showing that the SEm does not impact on the property of the PML while the DGm makes it unstable, even in isotropic medium where the continuous PML system has been proven to be stable [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]. This result is clearly in favor of coupling SEm with DGm, even in cases where the propagation domain justifies the use of triangles or tetrahedra. We go further by addressing the case of a heterogeneous medium.

Numerical results on a heterogeneous domain

Here we assess the case of the heterogeneous domain presented in Chapter 5 in Figure 5.13. We know that the PML formulation is unstable in anisotropic media. Hence we restrict ourselves to isotropic wave equation. We remind the equations we aim to solve on this domain: λ and µ are the Lamé parameters which have been defined in Chapter 2. The physical parameters V p , V s and ρ are the same as the ones given in Table 5.7 in Chapter 5. The domain is discretized using 55368 unstructured triangles and nearly 200000 squares. The high number of squares is due to the choice of surrounding the computational domain with squares in order to apply PML on every side and we can not reduce the number of cells by changing the space step because we do not have implemented h-adaptivity. The considered mesh is depicted on Figure 6.13. We apply the DG-SEm with DGm and SEm Figure 6.13: Heterogeneous case meshed and surrounded by quadrangles for PML(The PML layers are between the red lines)

     ρ∂ t v f = -∇p 1 ρc 2 0 ∂ t p = -∇ • v f (6.36a) (6.36b)
both at order 3. The time scheme used is a Runge-Kutta 4 scheme with a time-step ∆t equals to 1e -5. The final time of simulation is T = 2s. For the PML layers, we use the C-PML formulation with a damping profile d i and α i , ∀i = x, z defined with formula (6.34) and (6.35). For these simulations, we use the following parameters: L = 300m in every direction, N = 2, m = 1, V max p = 5334m.s -1 , R c = 0.0000001% . For these simulation we use the same second order Ricker source as for the homogeneous cases with f peak = 10Hz and t peak = 0.12s. It is placed in the water layer at point (1500,2650) Figure 6.14 presents the obtained snapshots and we observe that the PML absorbs the waves correctly in every direction and without any instability. This heterogeneous case shows thus that as in the previous case, the SEm applied in the PML contributes to stabilize the computations. We have performed the computations during a long time and we do not observe any spurious solution inside the layer.

Numerical results on a 3D elastic domain

To end this section on PML, we address the case of a three dimensional elastic domain with a PML layer placed on the top of the computational domain.

We consider a 1000x1000x1400 m elastic domain. In this domain we have a Pwavespeed equal to 2000m.s -1 and a S-wavespeed equal to 1000m.s -1 . The density ρ is equal to 1kg.m -3 everywhere in the domain. This domain is discretized using a hybrid mesh composed of 95511 unstructured tetrahedra and 50000 hexahedra placed at the top of the domain. For the space discretization, we use the DG-SEm with DGm and SEm of order two. The time integration is done using a RK4 scheme with a time step ∆t equals to 1e -5. The PML is placed at the top of the domain, where we apply SEm with cuboids and it is 300m thick. We use the C-PML formalism with the following parameters to define the damping d z and the polynomial α z : L = 300m, N = 2, m = 1, V max p = 2000m.s -1 , R c = 0.0000001% and C = 1. We use the following initial conditions: v s (x, 0) = 0, σ(x, 0) = 0 and the boundary condition σn = 0 on the external boundary of the computational domain, including the end of the PML.

The propagation is generated by a second order Ricker point source define as: S = (1 -2π 2 f 2 peak ((tt peak ) 2 )e -π 2 f 2 peak t 2 with a central frequency f peak equal to 20Hz and a delay t peak equal to 0.06s. This source is placed at point (500,500,950). The simulation lasts 0.8s.

Snapshots of the numerical waves are presented in Figure 6.15. The first slice at the top of each picture corresponds to the interface between cuboidal mesh and tetrahedral mesh. We see that as in the two dimensional case, the computation is stable and the PML absorbs correctly the waves.

Conclusion

Besides illustrating the implementation of PML inside the DG-SEm, this chapter provides an interesting application of the DG-SEm for stabilizing the PML. Indeed, even in isotropic media, the PML turns out to be unstable when full DGm is applied. We aim at pursuing in this direction by performing more numerical experiments.

Chapter 7

Conclusion

In the field of seismic imaging, computational domains are complex and very large. Hence reducing the computational costs for solving wave problems is an important issue that requires the evolution of existing numerical methods or the creation of new methods. In this thesis, we have proposed to couple two numerical methods that have proven their efficiency for solving wave equations: the Spectral Element Method (SEm) and the Discontinuous Galerkin method (DGm). This coupling has been carried out on hybrid meshes composed of structured quadrangular (cuboid) cells and unstructured triangular (or tetrahedral) cells. The use of hybrid meshes is justifies by the fact that SEm is easier to implement on quadrangle (or cuboid grids) and we wanted to keep the geometrical flexibility of unstructured triangle (or tetrahedral) meshes in which DGm had been developed previously by the team Magique-3D. To reach our objective, we have first developed DGm in quadrangle meshes. By this way, we aimed at verifying that when applied in the same structured mesh, SEm performs faster than DGm. We actually observed this property and it is worth mentioning that if DGm performs slower than SEm, it always provides a better relative error. Then we have proposed a coupled formulation combining SEm and DGm. For that purpose, we have first considered the area in which SEm is applied as a macro-element for DGm. This approach required considering a mixed-primal formulation of DGm to end up with a uniform expression for SEm and DGm with numerical fluxes acting as coupling terms at the interface between "SE" computational area and "DG" computational area.

It is worth noting that the coupling is achieved by assuming that the numerical interface is flat. The formulation has been proven to be stable thanks to an energy analysis leading to a stability result. The DG-SEm has been validated in various situations including the solution of the 3D elasto-acoustic wave equation. 2D numerical experiments clearly show that the DG-SE is faster than the full DGm. We acknowledge that more 3D experiments should be performed to get the same conclusions in 3D. In particular, it would be interesting to compare the DG-SEm in 3D with full DGm and full SEm. We have also investigated the introduction of domain truncation conditions into the DG-SEm by considering a first-order ABC and ADE-PMLs. The ABC has been implemented in the formulation as a penalization term and we have performed simulations which turn out to be stable in time but some reflections are generated by the external boundary. Regarding the implementation of PMLs, we have observed that the application of SEm inside the PML contributes to stabilize the computation while DGm introduces important numerical pollution turning quickly into blow up of the numerical solution. Some more investigations are on going to complete the works that have been performed in this PhD thesis. Our next objectives are:

• Go further by performing 3D simulations with a comparison with full DGm and possibly SEm. This is a reasonable perspective given that the full DGm is already operational and the full SEm is not far from being.

• Include local time-stepping to take advantage of the local characteristics of the mesh; At first, we planned to adapt the time-step to each numerical method, given that we knew that the CFL number of SEm is larger than the one of DGm. However, the numerical tests that we have performed in 2D for comparing full SEm and full DGm have shown that the SEm keeps outperforming the full DGm when applied with the same time-step. Hence we plan to improve the DG-SEm by adding local-time stepping in the DGm.

• Use the DG-SEm for solving 2D direct elasto-acoustic problems in view of Full Waveform Inversion (FWI) in the time domain. We have seen that DG-SEm is faster than the full DGm which is currently used in the team in a FWI algorithm that is currently under development in the team Magique-3D. The use of a DG-SE propagator should decrease the computational time of the method and thus widen the field of applications of the inversion method. This task should be achieved in collaboration with Total in the framework of the Total-Inria agreement DIP.

• Extend the method to the use of general hybrid meshes. The method has been developed for particular meshes composed on structured quadrangles or cuboids and we have assumed that in the vicinity of the numerical interface, a quadrangle is the neighbour of only one triangle and a cuboid is the neighbour of not more than two tetrahedra. For that purpose, we aim at using the recent developments in [START_REF] Reberol | Maillages hex-dominants: génération, simulation et évaluation[END_REF] which permit the construction of hybrid meshes including pyramids. Mots-clés Galerkine Discontinue, Elements Spectraux, couplage, maillage hybride Abstract To solve wave equations in heterogeneous media with finite elements with a reasonable numerical cost, we couple the Discontinuous Galerkin method (DGm) with Spectral Elements method (SEm). We use hybrid meshes composed of tetrahedra and structured hexahedra. The coupling is carried out starting from a mixed-primal DG formulation applied on a hybrid mesh composed of a hexahedral macro-element and a sub-mesh composed of tetrahedra. The SEm is applied in the macro-element paved with structured hexahedrons and the coupling is ensured by the DGm numerical fluxes applied on the internal faces of the macro-element common with the tetrahedral mesh. The stability of the coupled method is demonstrated when time integration is performed with a Leap-Frog scheme. The performance of the coupled method is studied numerically and it is shown that the coupling reduces numerical costs while keeping a high level of accuracy. It is also shown that the coupled formulation can stabilize the DGm applied in areas that include Perfectly Matched Layers. Keywords Discontinuous Galerkin, Spectral Elements, coupling, hybrid meshes

Résumé

1. 1

 1 Typical data acquisition for seismic marine imaging . . . . . . . . . . . . . . 2.1 Illustration of P-wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Illustration of horizontal shear (SH) wave . . . . . . . . . . . . . . . . . . . 2.3 Illustration of vertical shear (SV) wave . . . . . . . . . . . . . . . . . . . . . 2.4 Scheme for isotropic medium (left), Transverse Isotropy medium with one plan and one symmetry axis (middle), and orthotropic medium with two symmetric plans (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 TTI tilt angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Typical data acquisition for seismic marine imaging . . . . . . . . . . . . . . 2.7 Ricker wave for f peak = 10Hz and t peak = 0.2 . . . . . . . . . . . . . . . . . 3.1 Numerical solutions on a 1D mesh (continuous and discontinuous) . . . . . 3.2 h-adaptivity on a 2D mesh in P 1 . . . . . . . . . . . . . . . . . . . . . . . . 3.3 p-adaptivity on a 2D mesh with P 1 and P 2 . . . . . . . . . . . . . . . . . . . 3.4 Passage to the reference element K in 2D . . . . . . . . . . . . . . . . . . . 3.5 Domain of computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Speed modulus of the P -wave propagation in time using DGm . . . . . . . 3.7 P-wave propagation in a mesh with 22500 quadrilaterals using DGm . . . . 3.8 Seismograms of Gar6more solution (in blue) and numerical solution using DGm (in red) applied in an acoustic domain discretized with 10000 quadrangles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Seismograms of Gar6more solution (in blue) and numerical solution using DGm (in red) applied in an acoustic domain discretized with 22500 quadrangle cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Speed modulus of a wave propagating in an elastic domain. Simulation is done with DGm used on a mesh composed of 10000 cells. . . . . . . . . . . 3.11 Speed modulus of a wave propagating in an elastic domain. Simulation is done with DGm used on a mesh composed of 22500 quadrangles. . . . . . . 3.12 Seismograms of Gar6more solution (in blue) and numerical solution using DGm (in red) applied in an elastic domain discretized with 10000 quadrangle cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 8

 8 Comparison between SEm of order 5 and DGm of order 3 . . . . . . . . . . 3.9 Penalized DGm of order 3 and SEm of order 3, acoustic domain, 10000 quadrangles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Same CFL in acoustic domain . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Same CFL in elastic domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 DGm and SEm of order 3, Leap Frog scheme, 10000 squared cells, bi-layered acoustic domain, staircase interface . . . . . . . . . . . . . . . . . . . . . . . 3.13 DGm and SEm of order 3 applied at equal CFL, Leap Frog scheme, 10000 squared cells, bi-layered acoustic domain, staircase interface . . . . . . . . . 5.1 Physical parameters used of 2D homogeneous domains . . . . . . . . . . . . 5.2 DGm vs DG-SEm in an acoustic domain. . . . . . . . . . . . . . . . . . . . 5.3 DGm vs DG-SEm in the elastic case. . . . . . . . . . . . . . . . . . . . . . . 5.4 Physical parameters of the stratified domain . . . . . . . . . . . . . . . . . . 5.5 Heterogeneous domain, DGm of order 2 in triangles and 8 in squares, DG-SEm with DG at order 2 and SEm at order 8, RK4 scheme in time . . . . . 5.6 Heterogeneous domain, DGm of order 3 in triangles and 5 in squares, DG-SEm with DG at order 3 and SEm at order 5, RK4 scheme in time . . . . . 5.7 Physical parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Heterogeneous domain, DGm of order 3, DG-SEm of order 3, RK scheme in time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Physical parameters used in the 3D elasto-acoustic domain . . . . . . . . . Chapter 1

Figure 2 .

 2 Figure 2.1 to 2.3 present the different types of waves which can appear in acoustic or elastic media and how they affect the medium of propagation.

Figure 2 . 1 :

 21 Figure 2.1: Illustration of P-wave

Figure 2

 2 Figure 2.3: Illustration of vertical shear (SV) wave

Figure 2

 2 Figure 2.4: Scheme for isotropic medium (left), Transverse Isotropy medium with one plan and one symmetry axis (middle), and orthotropic medium with two symmetric plans (right).

  Figure 2.5: TTI tilt angle

  Figure 3.1: Numerical solutions on a 1D mesh (continuous and discontinuous)

  Figure 3.5: Domain of computation

Figure 3

 3 Figure 3.6: Speed modulus of the P -wave propagation in time using DGm

Figure 3

 3 Figure 3.7: P-wave propagation in a mesh with 22500 quadrilaterals using DGm

Figure 3 .

 3 Figure 3.10: Speed modulus of a wave propagating in an elastic domain. Simulation is done with DGm used on a mesh composed of 10000 cells.

Figure 3 .Figure 3

 33 Figure 3.11: Speed modulus of a wave propagating in an elastic domain. Simulation is done with DGm used on a mesh composed of 22500 quadrangles.

Figure 3 .

 3 Figure 3.13: Seismograms of Gar6more solution (in blue) and numerical solution using DGm (in red) applied in an elastic domain discretized with 22500 quadrangle cells

  and 3.18 picture snapshots in the acoustic case and they are completed with a seismogram on Figures 3.19 and 3.20 representing the numerical and analytical solutions. The snapshots agree with the physics, the last one illustrating the use of a free surface condition on the external boundary. The seismogram recorded at point

Figure 3 .

 3 Figure3.17: Speed modulus propagation in acoustic domain using SEm in a mesh with 10000 quadrilaterals

Figure 3 .Figure 3

 33 Figure 3.18: Speed modulus propagation in acoustic domain using SEm in a mesh with 22500 quadrilaterals

Figure 3 .

 3 Figure 3.20: Seismograms of Gar6more solution (in blue) and numerical solution using SEm (in red) applied in an acoustic domain discretized with 22500 quadrangle cells

Figure 3 Figure 3 Figure 3

 333 Figure 3.21: Speed modulus propagation in elastic domain using SEm in a mesh with 10000 quadrilaterals

Figure 3

 3 Figure 3.24: Seismograms of Gar6more solution (in blue) and numerical solution using SEm (in red) applied in an elastic domain discretized with 22500 quadrangle cells

Figure 3 .

 3 Figure 3.25: P-wave propagation on a mesh with 10000 quadrilaterals using DGm without penalization terms

  Figure 3.26: Behavior of L 2 -error in terms of approximation degree: penalized DGm and SEm with optimal CFL number, acoustic domain, 10000 quadrangles.

Figure 3 .

 3 Figure 3.29: Speed modulus propagation, DGm of order 3, Leap Frog scheme, 10000 squared cells, bi-layered acoustic domain, staircase interface

Figure 3 .

 3 Figure 3.30: Speed modulus propagation, SEm of order 3, Leap Frog scheme, 10000 squared cells, bi-layered acoustic domain, staircase interface

Figure 3

 3 Figure 3.31: Speed modulus propagation, DGm of order 3, Leap Frog scheme, 10000 structured triangle cells, acoustic domain

  Figure 4.2 depicts an example of such a configuration. It is important to define the internal boundaries denoted by Γ s,int , with Γ s,int = Γ 2,int ∪ Γ 1,2 . Remark 1: By exploiting the DGm construction, we avoid the use of mortar elements and thus the use of Lagrangian pre-conditioner with a dedicated solver which could have 4.1. DG-SEM VARIATIONAL FORMULATION IN THE ELASTIC CASE Ω 1,h Ω 2,h

Figure 4 . 1 :

 41 Figure 4.1: Decomposition of Ω s,h into two domains: Ω 1,h paved with quadrangles and Ω 2,h paved with unstructured triangles

Figure 4

 4 Figure 4.2: Ω 1,h : quadrangle macro-element; Ω 2,h paved with triangles; in blue, Γ 1,2 supports the boundary

  Figure 5.1: Domain of computation

Figure 5 .

 5 Figure 5.3 and 5.4 present the numerical results respectively obtained with DG-SEm and full DGm in a full acoustic propagation domain.We observe that both simulations provide the same results. In particular, it is interesting to observe that the matching between quadrangles and triangles is transparent for the numerical waves and the change of numerical method does not affect their propagation in the sense that there are no

Figure 5

 5 Figure 5.3: DG-SEm based simulation: acoustic domain, order 3 in space, RK4 scheme in time.

Figure 5 Figure 5 Figure 5

 555 Figure 5.4: Full DGm based simulation: acoustic domain, order 3 in space, RK4 scheme in time.

Figure 5

 5 Figure 5.7: Elastic case with DG-SEm in space and RK4 scheme in time.

Figure 5 . 8 :

 58 Figure 5.8: Elastic case with full DGm in space and RK4 scheme in time.

  Figure 5.9: Stratified domain with a water layer at the top

Figure 5 .

 5 Figure 5.10: Stratified domain meshed with unstructured triangles and big squares

Figure 5 .

 5 Figure 5.11: Heterogeneous domain, DGm of order 2 and SEm of order 8, RK4 in time.

Figure 5 .

 5 Figure 5.12: Heterogeneous domain, DGm of order 2 in the triangles and DGm of order 8 in the squares, RK4 in time.

Figure 5

 5 Figure 5.13: Topography of the salt dome Figure 5.16 presents the numerical simulation using the DG-SEm and Figure5.17 presents the result using full DGm. We see that we get the same results for both simulations. In particular, the combination of the two different meshes does not generate any numerical artefact. Figures 5.19 and 5.18 provide a superimposition of DG-SEm and DGm solutions with the reference solution which shows that the numerical solutions all conform correctly at the receiver. This is confirmed by the L 2 -error values presented in Table5.8. The computational times are still in favor of the DG-SEm coupling with a decrease by a factor of about 2. As in most of the experiments we have carried out, we still observed

  Figure5.16 presents the numerical simulation using the DG-SEm and Figure5.17 presents the result using full DGm. We see that we get the same results for both simulations. In particular, the combination of the two different meshes does not generate any numerical artefact. Figures 5.19 and 5.18 provide a superimposition of DG-SEm and DGm solutions with the reference solution which shows that the numerical solutions all conform correctly at the receiver. This is confirmed by the L 2 -error values presented in Table5.8. The computational times are still in favor of the DG-SEm coupling with a decrease by a factor of about 2. As in most of the experiments we have carried out, we still observed

Figure 5 .

 5 Figure 5.17: Heterogeneous domain, full DGm of order 3, RK4 in time.

Figure 5 .

 5 Figure 5.18: Seismograms of the reference solution (in blue) and numerical solution using DGm (in red).

VFigure 5

 5 Figure 5.19: Seismograms of the reference solution (in blue) and numerical solution using DG-SEm coupling (in red).

Figure 5 .

 5 Figure 5.20: Boundary configuration between an hexahedron and two tetraedra

Figure 5 .

 5 Figure 5.21: Cross section of the three dimensional mesh

Figure 5 .

 5 Figure 5.22: Wave propagation in a fluid-solid 3D domain using DG-SEm.

Figure 6 . 3 :

 63 Figure 6.3: Wave propagation in a fluid domain using DG-SEm with ABC. The three last snapshots are scaled up to see reflections at grazing incidence

Figure 6 . 4 :

 64 Figure 6.4: Wave propagation in a fluid domain using DGm with ABC. The three last snapshots are scaled up to see reflections at grazing incidence.

Figure 6 . 5 :

 65 Figure 6.5: Wave propagation in an isotropic elastic domain using DG-SEm with ABC. The three last snapshots are scaled up to see reflections at grazing incidence

Figure 6 . 6 :

 66 Figure 6.6: Wave propagation in an isotropic elastic domain using DGm with ABC. The three last snapshots are scaled up to see reflections at grazing incidence

Figure 6 . 7 :

 67 Figure 6.7: Heterogeneous domain, DG-SEm of order 3, ABC on each side of the external boundary.

(

  dx+αx)τ dτ = a x e -(dx+αx)m∆t anda x = d x d x + α x (b x -1), b x = e -dx+αx ∆t. (6.30)

Figure 6 . 9 :

 69 Figure 6.9: Speed modulus of an elastic wave propagation done using DGm and PML formulation.

Figure 6 .

 6 Figure 6.10: Speed modulus of an elastic wave propagation done using DGm and CPML formulation

Figure 6 .

 6 Figure 6.11: Speed modulus of an elastic wave propagation done using DG-SEm and PML formulation

Figure 6 .

 6 Figure 6.12: Speed modulus of an elastic wave propagation done using DG-SEm and CPML formulation

Figure 6 .

 6 Figure 6.14: Wave propagation on a domain representing a salt dome with PML all around the domain.

Figure 6

 6 Figure 6.15: Speed modulus, DGm and SEm at order 2, 3D elastic domain, PML layer at the top

  3.13 Seismograms of Gar6more solution (in blue) and numerical solution using DGm (in red) applied in an elastic domain discretized with 22500 quadrangle cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Position of the Gauss-Lobatto points in Q 4 . . . . . . . . . . . . . . . . . . 3.15 Passage to the reference element K in 2D . . . . . . . . . . . . . . . . . . . 3.16 Elements with local Gauss-Lobatto points(in red) and degrees of freedom on the mesh (in black) in Q 3 . Note that the local points inside the elements coincide with the global one. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Speed modulus propagation in acoustic domain using SEm in a mesh with 10000 quadrilaterals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Speed modulus propagation in acoustic domain using SEm in a mesh with 22500 quadrilaterals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Seismograms of Gar6more solution (in blue) and numerical solution using SEm (in red) applied in an acoustic domain discretized with 10000 quadrangle cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .27 Behavior of L 2 -error in terms of CFL number: penalized DGm at order 3, SEm at order 3 and SEm at order 5, acoustic domain, 10000 quadrangles. . 3.28 Bi-layered acoustic domain with a stair-like interface. . . . . . . . . . . . . 3.29 Speed modulus propagation, DGm of order 3, Leap Frog scheme, 10000 squared cells, bi-layered acoustic domain, staircase interface . . . . . . . . . 3.30 Speed modulus propagation, SEm of order 3, Leap Frog scheme, 10000 squared cells, bi-layered acoustic domain, staircase interface . . . . . . . . . Wave propagation in an isotropic elastic domain using DG-SEm with ABC. The three last snapshots are scaled up to see reflections at grazing incidence 6.6 Wave propagation in an isotropic elastic domain using DGm with ABC. The three last snapshots are scaled up to see reflections at grazing incidence 6.7 Heterogeneous domain, DG-SEm of order 3, ABC on each side of the external boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Position of the point source and the PML layer (in blue), homogeneous elastic domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Speed modulus of an elastic wave propagation done using DGm and PML formulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 Speed modulus of an elastic wave propagation done using DGm and CPML formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.11 Speed modulus of an elastic wave propagation done using DG-SEm and PML formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Correspondance table using Voigt notation . . . . . . . . . . . . . . . . . . 3.1 General Butcher table . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	List of Tables
	2.1

3.20 Seismograms of Gar6more solution (in blue) and numerical solution using SEm (in red) applied in an acoustic domain discretized with 22500 quadrangle cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Speed modulus propagation in elastic domain using SEm in a mesh with 10000 quadrilaterals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.22 Speed modulus propagation in elastic domain using SEm in a mesh with 22500 quadrilaterals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.23 Seismograms of Gar6more solution (in blue) and numerical solution using SEm (in red) applied in an elastic domain discretized with 10000 quadrangle cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.24 Seismograms of Gar6more solution (in blue) and numerical solution using SEm (in red) applied in an elastic domain discretized with 22500 quadrangle cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.25 P-wave propagation on a mesh with 10000 quadrilaterals using DGm without penalization terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.26 Behavior of L 2 -error in terms of approximation degree: penalized DGm and SEm with optimal CFL number, acoustic domain, 10000 quadrangles. . . . 33.31 Speed modulus propagation, DGm of order 3, Leap Frog scheme, 10000 structured triangle cells, acoustic domain . . . . . . . . . . . . . . . . . . . 4.1 Decomposition of Ω s,h into two domains: Ω 1,h paved with quadrangles and Ω 2,h paved with unstructured triangles . . . . . . . . . . . . . . . . . . . . . 6.5 6.12 Speed modulus of an elastic wave propagation done using DG-SEm and CPML formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.13 Heterogeneous case meshed and surrounded by quadrangles for PML(The PML layers are between the red lines) . . . . . . . . . . . . . . . . . . . . . 6.14 Wave propagation on a domain representing a salt dome with PML all around the domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.15 Speed modulus, DGm and SEm at order 2, 3D elastic domain, PML layer at the top . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Butcher table in the case of explicit RK scheme . . . . . . . . . . . . . . . . 3.3 Butcher table for RK2 scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Butcher table for RK4 scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Acoustic case, DG scheme without penalization, DGm and SEM of order 3, 10000 quadrangle cells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 CFL and CPU time comparison: acoustic domain, penalized DGm scheme. 3.7 CFL and CPU time comparison: elastic domain, penalized DGm scheme. .

  Table 3.6 and thus increase the number of time iterations. It represents a clear difference of CPU time that is clearly in favour of SEm.

		CFL	L2-error CPU-time Nb of time steps Number of cells
	DG	1e-3	1.259e-3	7.93	1000	10000
		6.66e-4 5.845e-4	32.98	1502	22500
	SEM 2.45e-3 7.903e-3	0.80	409	10000
		1.63e-3 3.608e-3	3.54	613	22500
	Table 3.6: CFL and CPU time comparison: acoustic domain, penalized DGm scheme.

Table 3 .

 3 7: CFL and CPU time comparison: elastic domain, penalized DGm scheme.

	CFL	L2-error CPU-time (s) Nb of time steps Number of cells
	DG 2.06e-3 1.071e-2	11.10	971	10000
	1.37e-3 4.731e-3	38.62	1456	22500
	SEM 6.09e-3 9.151e-2	1.19	356	10000
	4.06e-3 4.130e-2	4.18	534	22500

Table 3

 3 Behavior of L 2 -error in terms of CFL number: penalized DGm at order 3, SEm at order 3 and SEm at order 5, acoustic domain, 10000 quadrangles.

	.11.
	L2-error CPU-time Nb of cells
	DG 1.259e-3	7.93	10000
	4.731e-3	32.98	22500
	SEM 1.605e-3	2.12	10000
	4.756e-3	8.67	22500

.10: Same CFL in acoustic domain We can observe that in each case, SEm is much more efficient than DG in terms of computational time. Moreover, we can observe that SEm outperforms less DGm in the

Table 3 .

 3 11: Same CFL in elastic domain elastic case than in the acoustic case. Actually, if SEm is four times faster than DGm in a fluid domain, it is only twice faster in the solid domain. This is clearly in favour of the use of SEm in the fluid domain coupled with DGm in the solid area which may contain topography effects easier to handle with triangle cells. Hence, in the following of this work, the SEm will be used for solving the acoustic wave equation.

Table 3

 3

			L2-error CPU-time
	DG	1e-3	2.204e-3	16.04
	SEM 2.45e-3 1.280e-2	1.70
	Table 3.12: DGm and SEm of order 3, Leap Frog scheme, 10000 squared cells, bi-layered
	acoustic domain, staircase interface
		CFL L2-error CPU-time
	DG	1e-3 2.204e-3	16.04
	SEM 1e-3 1.458e-3	4.30

.13: DGm and SEm of order 3 applied at equal CFL, Leap Frog scheme, 10000 squared cells, bi-layered acoustic domain, staircase interface

  It is worth noting that in this

		V p (m/s) Vs (m/s) ρ (kg/m 3 )
	Area 1	1500	0	1
	Area 2	1800	1040	1.6
	Area 3	2100	1212	2
	Area 4	2600	1501	2.6
	Area 5	3100	1789	3.2
	Area 6	4000	2310	4
	Table 5.4: Physical parameters of the stratified domain

  Pour résoudre des équations d'ondes posées dans des milieux hétérogènes avec des éléments finis et un coût numérique raisonnable, nous couplons la méthode Discontinue de Galerkine (DGm) avec des éléments finis spectraux (SEm). Nous utilisons des maillages hybrides composés de tétraèdres et d'hexaèdres structurés. Le couplage est réalisé en partant d'une formulation DG mixte primale posée dans un maillage hybride composé d'un macro-élément hexaédrique et d'un sous-maillage composé de tétraèdres. La SEm est appliquée dans le macro-élément découpé en hexaèdres structurés et le couplage est assuré par les flux numériques de la DGm appliqués sur les faces internes du macroélément communes avec le maillage tétraédrique. La stabilité de la méthode couplée est démontrée quand l'intégration en temps est effectuée avec un schéma Saute-Mouton. Les performances de la méthode couplée sont étudiées numériquement et on montre que le couplage permet de réduire les coûts numériques avec un très bon niveau de précision. On montre aussi que la formulation couplée peut stabiliser la méthode DG appliquée dans des domaines incluant des couches parfaitement adaptées.

Γ s,h [[(Cξ)v s ]]. (3.19)Now we separate the terms on the boundary between external and internal boundaries:             Γ s,h [[σw s ]] = Γ s,h,out [[σw s ]] + Γ s,h,int [[σw s ]], Γ s,h [[(Cξ)v s ]] = Γ s,h,out [[(Cξ)v s ]] + Γ s,h,int [[(Cξ)v s ]].

Remerciements

that the relative error is lower for DGm than for SEm with a reduction by a factor of about 40. Table 5.8: Heterogeneous domain, DGm of order 3, DG-SEm of order 3, RK scheme in time Now that we have seen the efficiency of the coupling in the 2D case, we illustrate how DG-SEm performs for solving the 3D elasto-acoustic system.

3D elasto-acoustic wave equation

The extension to 3D problems raises some difficulties regarding the construction of hybrid meshes. In this work, we restrict the study to the case of meshes composed of cuboids and tetrahedra and a cuboid can only have two tetrahedra as neighbours (see Figure 5.20 for an example).

We present a simulation done in a 3D elasto-acoustic domain. The domain is of size 1000x1000x1400. The water layer is on the top of the computational domain and is 400 meters thick. The physical interface between water and subsurface coincides with the numerical interface. The physical parameters associated with this case are given in Table 5.9.

By this way, it is possible to use time scheme of order higher than 2 (like Runge-Kutta 4 for example). Moreover, it has been proved in [START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML)[END_REF] that the ADE-PML formulation is equivalent to the C-PML formulation when using a second-order time scheme. Here we describe the method by restricting ourselves to the first equation of (6.21), the others variables are handled in exactly the same way. This equation reads:

Then according to (6.24), we have:

Now we define the auxiliary memory variable noted Ψ p,x : (6.32) Using this definition, equation (6.31) becomes:

and Ψ px satisfies the following equation :

Then by going back to the time domain, we obtain:

To construct the ADE-PML formulation, we then have to introduce an auxiliary variable Ψ z associated with v f,z and in the same way, the pressure is associated with two auxiliary unknowns. Finally we have the following ADE-PML system:

We end up with a system of 7 variables which has the following discrete form:

where the matrices M v f,h , M p h and R p h , R v f,h,x and R v f,h,z are the mass and stiffness matrices defined in chapter 3.