Marie Bonnasse-Gahot

Henri Calandra

Julien Diaz

Stéphane Lanteri

Discontinuous Galerkin methods for the simulation of the propagation of the elastic wave equations in the frequency domain

Keywords: seismic imaging, forward problem, frequency domain, elastic waves, centered flux DG method, upwind flux DG method, discontinuous Galerkin method imagerie sismique, résolution du problème direct, domain fréquentiel, ondes élastiques, méthodes de Galerkin discontinues, méthode GD à flux centrés, méthode

In this report, we study two nodal discontinuous Galerkin (DG) methods, the centered flux DG method and the upwind flux DG method, for the resolution of the 2D elastic waves equations in frequency domain, so called Helmholtz equations. We give the formulation of both methods and we compare the obtained results.

Introduction

The numerical solution of wave equations (acoustics, elastodynamics or electromagnetics) is increasingly used in many areas such as civil engineering, mechanical engineering, aerospace engineering, geophysics, medicine, biology or telecommunications. Applications are then various: noise reduction, radar and antenna design, satellites and waveguides, detection of hidden targets, medical imaging, seismic imaging, earthquakes dynamics, etc.

The scientific context of the present study is seismic imaging, which aims at recovering the structure of the earth. A seismic imaging method using a wave equation model is called a seismicreflection technique. It is the most used imaging method in the petroleum industry because it is the one that yields the most accurate images of the subsurface in view of analyzing them to find hydrocarbons.

The principle of a seismic acquisition is simple: sources, generally placed on the top of the subsurface, emit waves inside the earth; then receivers, which can be placed at the top or in the depth, record the reflected waves (arrival time and amplitude). The amplitudes of reflected fields allow to recover the material characteristics constituting the ground, and the arrival times allow to establish where the reflectors are. Reflectors are the interfaces between two different media.

Seismic imaging campaigns can be done in the sea or in the ground and with its results, one can construct a card that represents the variations of the velocity in the medium. The latter is referred to as the velocity model. Its quality depends on the number of sources used. It is the reason why in seismic imaging the number of sources used is usually large (about 1000). Consequently, the efficiency of the whole procedure is directly related to the efficiency of the numerical method used to solve the wave equations.

There exists several seismic imaging methods: the Reverse Time Migration (RTM) and the Full Wave Inversion (FWI) are the most widely used. Both methods are based on the resolution of wave equations.

The RTM uses the fact that the wave equations are reversible. From the sources, one propagates waves to obtain the propagated field. Then one makes use of the recorded data of the reflected waves as initial data, and propagates them in order to obtain the retro-propagated field. Finally the image of the subsurface is obtained by correlating the propagated fields and retro-propagated fields: at each point where there is a correlation, a reflector is deduced.

The FWI is an iterative procedure solving 2N harmonic wave equations at each iteration of the algorithm if N sources are used. It defines an inversion process (see [START_REF] Kern | Problèmes inverses: aspects numériques[END_REF] for more details about inversion process), and thus it is composed of two major steps: first, one solves the forward problem (i.e. the propagation of the wave equations) which models the phenomenon for each source; then, one computes residuals thanks to a comparison between the recorded data and the numerical data using, e.g., a least-squares method. These residuals are used to update the velocity and the process is repeated until there is accordance between recorded data and numerical data.

Seismic imaging can be performed in the time-domain or in the frequency-domain. Timedomain approaches do not require an important computational cost, but the implementation of the imaging condition is more complicated than in the frequency-domain case. Furthermore, in frequency-domain approaches, it is not necessary to store the solution at each time step of the forward simulation. This is interesting because seismic imaging involves very large problems with a lot of data. Memory must be used with attention. In the frequency-domain case, the main drawback then lies in the need to solve a large linear system of equations which represents a challenging task when considering realistic 3D elastic media, despite the recent advances on high performance computing.

In this context, we are interesting in the modeling part, i.e. the resolution of the forward problem of the FWI, assuming a time-harmonic regime, leading to the so-called Helmholtz equations. The main objective of our work is to propose and develop a new finite element type solver characterized by a reduced-size discrete operator (as compared to existing such solvers) without hampering the accuracy of the numerical solution.

A wide variety of approximation methods of the wave equations is currently available. Moreover, with the progress of high performance computing, it is now possible to model elastic domain. The FWI is performed with methods belonging to the family of finite difference (FD) method or finite element (FE) methods. FD methods use regular grid and allow to obtain easier systems to solve (as compared to FE methods). But they may lack of sufficient accuracy when one considers a highly heterogeneous medium or when the underlying mesh is too coarse. Moreover with FD methods, one does not take into account the irregular topography of the subsurface because the structured mesh can not correctly approach interfaces. Indeed, with Cartesian meshes, it is difficult to handle a steep subsurface.

By contrast, with FE methods, one can use unstructured meshes to discretize accurately complex domains. One inherits from a greater flexibility in the construction of the mesh. However, FE methods require more memory space than FD methods. Two kinds of FE methods seem to be adapted to wave propagation simulation: spectral elements (SE) methods and discontinuous Galerkin (DG) methods (see [START_REF] Virieux | Modelling seismic wave propagation for geophysical imaging[END_REF] for a comparison between DG methods and finite differences methods applied to wave equations for seismic applications).

Inria SE methods use high order functions and need less memory space than classical FE methods without hampering the numerical convergence order. However, SE methods are formulated on quadrangular (2D case)/hexahedral (3D case) meshes which may not be ideally adapted to the discretization of complex geometries in the 3D case. DG methods and FE methods mainly differ on basis functions: DG basis functions are only piecewise continuous. Moreover in addition to the fact that they are formulated on unstructured triangular meshes, they are more suitable than continuous Galerkin (CG) methods to deal with hp-adaptivity (interpolation degree p or mesh step h can change from element to another), providing a greater flexibility in the mesh construction and the discretization of the different regions of the subsurface. For time-domain problems, DG methods provide explicit representation of the solution because the mass matrix is block diagonal. In addition, they are nicely adapted to high performance computing.

With the goal to develop more accurate, more efficient and cheaper forward solvers, we are interested here in solving the elastic Helmholtz equations with DG methods. In his doctoral thesis, R. Brossier [START_REF] Brossier | Imagerie sismique à deux dimensions des milieux visco-élastiques par inversion des formes d'ondes: développements méthodologiques et applications[END_REF] applied the centered flux DG method to the elastic Helmholtz equations under a pseudo-conservative form. In this study we study numerically a centered flux DG method which is an adaptation of the method studied by S. Delcourte et al. in [START_REF] Delcourte | A high order discontinuous Galerkin method for the seismic wave propagation[END_REF] in the time-domain case, and an upwind flux DG method formulated as in [START_REF] Lanteri | High order discontinuous Galerkin method for the solution of 2D time-harmonic Maxwell's equations[END_REF] for the time-harmonic Maxwell equations. We note that Käser and Dumbser 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] have also considered an upwind flux DG formulation for time-domain elastic wave equations.

This report describes a preliminary study in view of the development a new class of DG methods, the hybridizable DG method, which will be the subject of a second document where we will in particular compare results obtained with upwind DG method and the new HDG method.

2 Problem statement and notations

Isotropic elastodynamics system

We consider the first order formulation of the elastodynamics system assuming a timeharmonic regime. In the 3D case, for x = (x, y, z) ∈ Ω ⊂ R 3 , this systems reads as

iωρ(x)v(x) = ∇ • σ(x) + f (x) in Ω, iωσ(x) = C(x) (v(x)) in Ω, (2.1)
where i is the imaginary unit, ω the angular frequency, ρ(x) is the mass density and f (x) is the source term (volumetric forces). Moreover, v(x) = (v x (x), v y (x), v z (x) T is the velocity vector, is the stress tensor, with ij = 1 2 ∂v i ∂j + ∂v j ∂i , i, j = x, y, z, σ is the stain tensor with, in the isotropic case, σ ij = λδ ij Tr( ) + 2µ ij , i, j = x, y, z and C is a 3 × 3 × 3 × 3 symetric fourth order tensor, with elastic coefficients. Using Voigt's notation, C can be reduced to 6 × 6 matrix. Using

ij → α or kl → β =11 → 1 22 → 2 33 → 3 12 → 4 23 → 5 13 → 6
we have that in the isotropic case

C(x) =         λ(x) + 2µ(x) λ(x) λ(x) 0 0 0 λ(x) λ(x) + 2µ(x) λ(x) 0 0 0 λ(x) λ(x) λ(x) + 2µ(x) 0 0 0 0 0 0 µ(x) µ(x) 0 0 0 0 0 µ(x) µ(x) 0 0 0 µ(x) 0 µ(x)        
, where λ and µ are the Lamé's coefficients.

The boundary conditions are given by Thereafter, to simplify, we omit the spatial-dependency in the physical parameters ρ, λ and µ, Inria tensors C, σ and and vector v. We can develop the equations of system (2.1) and obtain (2.3)

σn = 0 on Γ a , (D n ) -W = 0 on Γ b , (2.2
                                                               iωV x =
In 2D, the y components do not exist and we only have a system of 5 equations with 5 unknowns

                             iωV x = 1 ρ ∂σ xx ∂x + ∂σ xz ∂z , iωV z = 1 ρ ∂σ xz ∂x + ∂σ zz ∂z , iωσ xx = (λ + 2µ) ∂V x ∂x + ∂V z ∂z , iωσ zz = (λ + 2µ) ∂V z ∂z + ∂V x ∂x , iωσ xz = µ ∂V x ∂z + ∂V z ∂x .
(2.4)

Vectorial form

We can rewrite system (2.3) under a vectorial form

iωW + A x ∂W ∂x + A y ∂W ∂y + A z ∂W ∂z = 0, (2.5)
with

W = (V x , V y , V z , σ xx , σ yy , σ zz , σ xy , σ yz , σ xz ) T , A x = -                    0 0 0 1 ρ 0 0 0 0 0 0 0 0 0 0 0 1 ρ 0 0
0 0 0 0 0 0 0 0 1 ρ λ + 2µ 0 0 0 0 0 0 0 0 λ 0 0 0 0 0 0 0 0 λ 0 0 0 0 0 0 0 0 0 µ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 µ 0 0 0 0 0 0

     0 0 0 0 0 0 1 ρ 0 0 0 0 0 1 ρ 0 0 0 0 0 0 0 0 0 0 1 ρ 0 λ 0 0 0 0 0 0 0 λ + 2µ 0 0 0 0 0 0 0 λ 0 0 0 0 0 0 µ 0 0 0 0 0 0 0 0 0 µ 0 0 0 0 0 0 0 0 0 0 0 0 0      0 0 0 0 0 0 0 0 1 ρ 0 0 0 0 0 0 0 1 ρ 0 0 0 0 0 0 1 ρ 0 0 0 0 0 λ 0 0 0 0 0 0 0 0 λ 0 0 0 0 0 0 0 0 λ + 2µ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 µ 0 0 0 0 0 0 0 µ 0 0 0 0 0 0 0 0                    .

Notations

We discretize the computational domain Ω in a collection of disjoint elements (tetrahedrons for the 3D case, triangles for the 2D case) T h . We denote by • F(K) the set of faces (edges in the 2D case, anyhow, we call it face) of an element K,

• F a face of an element K,

• F b the union of all boundary faces F b , i.e. F b = ∂K ∩ ∂Ω,

• F i the union of all interior faces F i i.e. F i = ∂K ∩ ∂K where K and K are neighbourgs,

• F h the set of all faces T h , i.e.

F h = F i ∪ F b ,
• n the outward unit norm vector to K, n the outward unit norm vector to K .

We define the jump [[•]] and the mean (average) {•} such as

• On a face F ∈ F i , the jump [[•]] of a vector u is [[u • n]] = u K • n + u K • n = u K • n -u K • n,
and the jump of a tensor σ is

[[σ • n]] = σ K n + σ K n = σ K n -σ K n. Inria • On a face F ∈ F b , the jump [[•]] of u is [[u • n]] = u K • n,
and for σ [[σ • n]] = σ K n.
• On a face F ∈ F i , the average {•} of a scalar variable u is

{u} = u K + u K 2 ,
and for a vector u

{u} = u K + u K 2 .
• On a face F ∈ F b , we have for u {u} = u K , and for u {u} = u K .

Let P p (D) denote the space of polynomial functions of degree at most p on domain D. For any element K of T h , Φ K = ϕ K i , i = 1, df K , is a basis of P p (K) with Lagrange's polynomial functions; df K is the number of degrees of freedom (ndof) with df K = (p + s)! s! with s the space dimension, and

ϕ K i ∈ P p , ϕ K i (a j ) = δ ij , 1 ≤ j ≤ df K ,
where δ ij is Kronecker symbol and a j coordinates of the j th K dof. We define the following finite elements spaces

V h = {v h ∈ L 2 (Ω)|∀K ∈ T h , v h | K ∈ P p (K)}, V h = {v h ∈ L 2 (Ω) 3 |∀K ∈ T h , v h | K ∈ (P p (K)) 3 }, Σ h = {σ h ∈ L 2 (Ω) 6 |∀K ∈ T h , σ h | K ∈ (P p (K)) 6 },
where L 2 (Ω) is the space of square integrable functions on the domain Ω.

3 Discontinuous Galerkin method

3.1 DG formulations

Principles

We decide to work with the vectorial equation (2.5) of the elastodynamics system in order to establish the DG formulations considered in this study. Moreover, we choose to make a classic approximation of coefficients ρ, λ and µ considering them as piecewise-constant, i.e. constant over an element. Matrices A x , A y and A z , depending only on these coefficients, are also piecewiseconstant.

The DG method seeks an approximate solution

W h = v h , σ h T in the space V h × Σ h satisfying for all K in T h K iωW h Φ + K A x ∂W h ∂x Φ + K A y ∂W h ∂y Φ + K A z ∂W h ∂z Φ = 0,
where Φ ∈ V h . An integration by parts yields

K iωW h Φ - K A x W h ∂Φ ∂x + ∂K A x W h Φn x - K A y W h ∂Φ ∂y + ∂K A y W h Φn y - K A z W h ∂Φ ∂z + ∂K A z W h Φn z = 0. (3.1)
Defining matrix D n such as

D n = n x A x + n y A y + n z A z = -                   0 0 0 n x ρ 0 0 n y ρ 0 n z ρ 0 0 0 0 n y ρ 0 n x ρ n z ρ 0 0 0 0 0 0 n z ρ 0 n y ρ nx ρ n x (λ + 2µ) n y λ n z λ 0 0 0 0 0 0 n x λ n y (λ + 2µ) n z λ 0 0 0 0 0 0 n x λ n y λ n z (λ + 2µ) 0 0 0 0 0 0 n y µ n x µ 0 0 0 0 0 0 0 0 n z µ n y µ 0 0 0 0 0 0 n z µ 0 n x µ 0 0 0 0 0 0                   = Dn, with D = (A x , A y , A z ) .
We rewrite (3.1) as

K iωW h Φ - K A x W h ∂Φ ∂x - K A y W h ∂Φ ∂y - K A z W h ∂Φ ∂z + ∂K (D n W h ) | ∂K Φ = 0. (3.2)
We note that the term (D n W h ) | ∂K will require a special treatment since W h is discontinuous at an element boundary and D n is constant on an element. This will be detailed later. By summing equations (3.2) over all elements, we obtain the following global equation

K∈T h K iωW h Φ - K∈T h K A x W h ∂Φ ∂x - K∈T h K A y W h ∂Φ ∂y - K∈T h K A z W h ∂Φ ∂z + F ∈Fi F D K n W K h + D K n W K h Φ + F ∈F b F D K n W K h Φ = 0. (3.3)
We choose to define the jump of D n W h through a face F such as

[[D n W h ]] = D K n W K h + D K n W K h for F ∈ F i , [[D n W h ]] = D K n W K h for F ∈ F b .
Finally the global DG formulation is given by

K∈T h K iωW h Φ - K∈T h K A x W h ∂Φ ∂x - K∈T h K A y W h ∂Φ ∂y - K∈T h K A z W h ∂Φ ∂z + F ∈F h F [[D n W h Φ]] = 0.
(3.4)

Inria

If the support of the basis function Φ is the element K, the global formulation (3.4) reduces to a local formulation on element K. The surface integral over K's boundary is written like the sum over all K faces

∂K D n W h Φ = F ∈F (K) F D n W h Φ.
Moreover as matrices A x , A y and A z are piecewise-constant, we obtain the following local equation

K iωW K h Φ K - K A K x W K h ∂Φ K ∂x - K A K y W K h ∂Φ K ∂y - K A K z W K h ∂Φ K ∂z + F ∈F (K) F (D n W h ) | F Φ K = 0. (3.5)
To compute the numerical trace (D n W h ) | F in equation (3.5), we consider two options: a centered scheme and an upwind scheme.

Centered flux DG scheme

To establish the centered flux DG formulation, we approximate

(D n W h ) | F on a face F by its average, i.e (D n W h ) | F = {D n W h } .
Using this definition, the local equation (3.5) becomes

K iωW K h Φ K - K A K x W K h ∂Φ K ∂x - K A K y W K h ∂Φ K ∂y - K A K z W K h ∂Φ K ∂z + F ∈∂K∩∂K F 1 2 D K n W K h + D K n W K h Φ K + F ∈∂K∩∂Ω F D K n W K h Φ K = 0. (3.6)
We verify that summing the centered flux local equation over all the elements yields the expected global DG equation

K∈T h K iωW K h Φ K - K∈T h K A K x W K h ∂Φ K ∂x - K∈T h K A K y W K h ∂Φ K ∂y - K∈T h K A K z W K h ∂Φ K ∂z + F ∈F h F {D n W h }[[Φ]] = 0.
Let us now recall the property linking the average and the jump

[[uv]] = {u}[[v]] + [[u]]{v},
and the fact that if W h is solution of the elastodynamics equations (2.1) then,[[W h ]] = 0. Thus, if we multiply W h by a constant matrix D n , we have [[D n W h ]] = 0, and we can deduce that

F h {D n W h }[[Φ]] = F ∈F h F {D n W h }[[Φ]] + F ∈F h F [[D n W h ]]{Φ} = F h [[D n W h Φ]].
Replacing the surface integral we recover the global equation (3.4).

In equation (3.6), we decompose each component of the vector W K h in the basis

Φ K = ϕ K i W K h = di j=1 W K j ϕ K j , (3.7)
where

W K j = v xj , v yj , v zj , σ xxj , σ xxj , σ yyj , σ zzj , σ xyj , σ yzj , σ xzj T .
Choosing the test function as

ϕ K l ∈ Φ K , l = 1, d i , we obtain the following local discretization iω di j=1 W K j K ϕ K j ϕ K l -A K x di j=1 W K j K ϕ K j ∂ϕ K l ∂x -A K y di j=1 W K j K ϕ K j ∂ϕ K l ∂y -A K z di j=1 W K j K ϕ K j ∂ϕ K l ∂z + F ∈∂K∩∂K 1 2   D K n di j=1 W K j F ϕ K j ϕ K l + D K n di j=1 W K j F ϕ K j ϕ K l   + F ∈∂K∩∂Ω F D K n W K h ϕ K l = 0. (3.8)
Finally we write the above equation in a matrix form

iωM K W K -A K x D K x W K -A K y D K y W K -A K z D K z W K + F ∈∂K∩∂K 1 2 D K n E K 1 W K + D K n E K 2 W K + F ∈∂K∩∂Ω F D K n W K h Φ K = 0, (3.9)
where

W K = W K 1 , ..., W K di T , M K (i, j) = K ϕ K j ϕ K i , D K u (i, j) = K ϕ K j ∂ϕ K i ∂u , with u = x, y, z, E K 1 (i, j) = F ϕ K j ϕ K i , E K 2 (i, j) = F ϕ K j ϕ K i .
(3.10)

Inria

We note that we have not developed the boundary integral

F ∈∂K∩∂Ω F D K n W K h Φ K .
We will detail the treatment of this term in section 3.1.4. By summing the local equation (3.9) over all elements, the global equation (3.4) can be written as

iω K∈T h M K W K - K∈T h A K x D K x W K - K∈T h A K y D K y W K - K∈T h A K z D K z W K + K∈T h F ∈∂K∩∂K 1 2 D K n E K 1 W K + D K n E K 2 W K + K∈T h F ∈∂K∩∂Ω F D K n W K h Φ K = 0.
(3.11)

Upwind flux DG scheme

Before developing the upwind flux DG formulation, we need to look at the some properties of the matrix D n . It can be shown that D n is a square matrix that can be diagonalized. The corresponding eigenvalues are

ξ 1 = -v p , ξ 2 = ξ 3 = -v s , ξ 4 = ξ 5 = ξ 6 = 0, ξ 7 = ξ 8 = v s , ξ 9 = v p , where v p = λ + 2µ ρ and v s = µ ρ
. The associated eigenvectors are gathered in 9 × 9 matrix R n , where kth column is the eigenvector associated of the eigenvalue

ξ k R n =               n x v p n z v s 0 0 0 0 • • • n y v p 0 n z v s 0 0 0 • • • n z v p -n x v s -n y v s 0 0 0 • • • λ + 2 µ n 2 x 2 µ n x n z 0 2 n y n z 0 0 • • • λ + 2 µ n 2 y 0 2 µ n y n z 0 2 n x n z 0 • • • 2 µ n 2 z + λ -2 µ n x n z -2 µ n y n z 0 0 2 n x n y • • • 2 µ n y n x µ n y n z µ n x n z -n x n z -n y n z n 2 z • • • 2 µ n y n z -µ n x n y µ n 2 z -n 2 y n 2 x -n y n x -n z n x • • • 2 µ n z n x µ n 2 z -n 2 x -µ n x n y -n x n y n y 2 -n z n y • • •               ,               • • • 0 -n z v s -n x v p • • • -n z v s 0 -n y v p • • • n y v s n x v s -n z v p • • • 0 2 µ n x n z λ + 2 µ n 2 x • • • 2 µ n y n z 0 λ + 2 µ n 2 y • • • -2 µ n y n z -2 µ n x n z 2 µ n 2 z + λ • • • µ n x n z µ n y n z 2 µ n y n x • • • µ n 2 z -n 2 y -µ n x n y 2 µ n y n z • • • -µ n x n y µ n 2 z -n 2 x 2 µ n z n x               . We define D + n such as D + n = R n Λ + (R n ) -1 ,
and D - n such as

D - n = R n Λ -(R n ) -1 ,
where Λ + is the diagonal matrix built with the positive eigenvalues of D n and Λ -the one including negative eigenvalues of D n ; R n is the matrix defined above. We remark that we have

D n = D + n + D - n and |D n | = D + n -D - n .
(3.12)

We remind the local equation (3.5) of the DG formulation

K iωW K h Φ K - K A K x W K h ∂Φ K ∂x - K A K y W K h ∂Φ K ∂y - K A K z W K h ∂Φ K ∂z + F ∈F (K) F (D n W h ) | F Φ K = 0.
For the upwind flux DG formulation, we define the numerical trace

(D n W h ) | F of D n W on a face F as (D n W h ) | F = D K n + W K h + D K n - W K h . (3.13)
By replacing in equation (3.5) the surface term by its approximation for an interior face, we obtain

K iωW K h Φ K - K A K x W K h ∂Φ K ∂x - K A K y W K h ∂Φ K ∂y - K A K z W K h ∂Φ K ∂z + F ∈∂K∩∂K F D K n + W K h + D K n - W K h Φ K + F ∈∂K∩∂Ω F D K n W K h Φ K = 0.
Using the decomposition (3.7) of W K in the basis Φ K and taking ϕ K l , l = 1, d i as test-function, we obtain the local equation for the upwind DG schema

iω di j=1 W K j K ϕ K j ϕ K l -A K x di j=1 W K j K ϕ K j ∂ϕ K l ∂x -A K y di j=1 W K j K ϕ K j ∂ϕ K l ∂y -A K z di j=1 W K j K ϕ K j ∂ϕ K l ∂z + F ∈∂K∩∂K   D K n + di j=1 W K j F ϕ K j ϕ K l + D K n - di j=1 W K j F ϕ K j ϕ K l   + F ∈∂K∩∂Ω F D K n W K h ϕ K l = 0.
(3.14)

Inria

Similarly to the centered scheme, we write the above equation in a matrix form

iωM K W K -A K x D K x W K -A K y D K y W K -A K z D K z W K + F ∈∂K∩∂K D K n + E K 1 W K + D K n - E K 2 W K + F ∈∂K∩∂Ω F D K n W K h Φ K = 0, (3.15)
where matrices M K , D K u , u = x, y, z, E K 1 and E K 2 are defined by (3.10). Finally the global equation for the upwind DG scheme in a matrix form is

iω K∈T h M K W K - K∈T h A K x D K x W K - K∈T h A K y D K y W K - K∈T h A K z D K z W K + K∈T h F ∈∂K∩∂K D K n + E K 1 W K + D K n - E K 2 W K + K∈T h F ∈∂K∩∂Ω F D n Wϕ = 0.
(3.16)

Boundary conditions

We are concerned here with the numerical treatment of the integral term over the boundary of the computational domain,

F ∈∂K∩∂Ω F D n Wϕ.
Two boundary conditions have been introduced in the problem statement in section 2.1

• Free surface condition to simulate a physical interface between the domain and the air

• Absorbing boundary condition to simulate an infinite domain

Absorbing boundary condition. The objective is to define an appropriate numerical treatment when the theoretically unbounded propagation domain is artificially truncated. More precisely, we would like impose some conditions on the corresponding boundary faces which allow the absorption of the waves when they reach the artificial boundary. There exist several possible strategies including Perfectly Matched Layers (PML, C-PML, etc.) absorbing conditions which are approximations of an exact transparent condition. The latter option has be adopted in this study. By using the plus-minus decomposition of the flux matrix (3.12), we can define a simple treatment which consists in considering that the inflow flux is zero for each face on the artificial boundary, i.e.

(D n ) -W = 0 on F b . (3.17) Then F b D n W h ϕ F b D K n + W K h ϕ D K n + di j=1 W K j F b ϕ K j ϕ.
Setting ϕ = ϕ K i we have

F b D K n + W K h ϕ K i = D K n + W K B K , where B K (i, j) = F b ϕ K j ϕ K i .
Free surface condition. On a free surface boundary we impose

σn = 0, (3.18)
and no particular condition is applied to the velocity components. We then have

F b D K n W K h ϕ = F b D K n LW K h ϕ, where L =              
1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

             
. Using ϕ = ϕ K i the numerical flux on free surface boundary is

F b D K n LW K h ϕ K i = D K n LW K B K .
Finally, the global discrete system writes iω

K∈T h M K W K - K∈T h A K x D K x W K - K∈T h A K y D K y W K - K∈T h A K z D K z W K + K∈T h F ∈∂K∩∂K H K n E K 1 W K + H K n E K 2 W K + K∈T h F ∈∂K∩∂Ω G n B K W K = 0, (3.19)
where

H K n = D K n + and H K n = D K n - in the case of the upwind numerical flux, H K n = 1 2 D K n and H K n = 1 2 D K n for the centered numerical scheme, and G n = D K n L or G n = D K n +
depending on the type of boundary conditions.

Two-dimensional numerical results

In this section, we provide some numerical results in 2D to assess the performances (accuracy and efficiency) of the centered and upwind DG schemes introduced in the previous section. These Inria schemes have been implemented in a Fortran 90 software. We use the MUMPS sparse direct solver for the resolution of the linear system (see [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF] for more details) resulting from the DG discretization scheme.

Numerical experiments are performed on a hardware system equipped with 2 quad-core Nehalem Intel Xeon X5550/2,66 GHz CPUs and 24Go (DDR3 1333 MHz) of RAM.

Two simple problems, the propagation of a plane wave in a homogeneous medium and the scattering of a plane wave by a disk, and one heterogeneous problem, the scattering of a plane wave by a solid disk, are considered. We remind that to propagate the waves we have to solve the elastic Helmholtz equations

iωρ(x)v(x) = ∇ • σ(x) + f (x)
in Ω,

iωσ(x) = C(x) (v(x)) in Ω, (4.1)

Plane wave propagation in an homogeneous medium

We first consider the simple test problem of the propagation of a plane wave in an homogeneous medium. The computational domain Ω is a 10000 m × 10000 m square. The physical properties of the medium are ρ = 1000 kg.m -3 and values of Lamé's coefficients λ and µ that are set to 8 M P a and 4 M P a respectively. These values imply a velocity v p of P -waves equal to 4000 m.s -1 and a velocity v s of S-waves equal to 2000 m.s -1 . On the boundaries we impose an absorbing condition with a plane wave incident field

U = ∇e i(kx cos θx+kz sin θz) =       V x0 V z0 σ xx0 σ zz0 σ xz0       e i(kx cos θx+kz sin θz) , where k = k 2 x + k 2 z = ω v p
is the wavenumber, k x = ω v p cos θ and k z = ω v p sin θ, and θ is the incidence angle. ω is the angular frequency, ω = 2πf where f is the frequency. If we choose arbitrarily V x0 , we can express the other components such as

                       V z0 = k x k z (λ + µ) ρω 2 -k 2 x µ -k 2 z λ + 2µ) V x0 σ xx0 = -1 ω (k x (λ + 2µ) V x0 + λk z V z0 ) σ zz0 = -1 ω (λk x V x0 + (λ + 2µ) k z V z0 ) σ xz0 = -µ ω (k z V x0 + k x V z0 ) .
In the simulations we choose θ equal to 0 and f = 2 Hz, then ω = 4π 12.56 Hz. We discretize the computational domain Ω into three unstructured meshes with respectively 3000, 10 000 and 45 000 elements. Two of these meshes are shown on figs. 4.1 and 4.2; their characteristics are given in table 4.1. Figure 4.9: Numerical solution, mesh M2, upwind flux DGFD-P 1 scheme, V x component. Fig. 4.14 shows the numerical convergence of the centered and upwind DGFD methods. These two graphs confirm that the centered DGFD scheme is less accurate than the upwind DGFD scheme and we observe that the centered DGFD scheme converges with order p whereas upwind DGFD scheme converges with order p + 1 i.e. with optimal rate. The computational performances of both methods for all simulations of this test problem are summarized in tables 4.2 for the number of non-zero terms in the global matrix and the memory used (in MB), 4.3 for the time required for the construction of the global matrix (in seconds) and for the resolution (in seconds) and 4.4-4.5 for the mean and relative errors in the V x and σ xx components. The upwind DG scheme requires more memory space than the centered DG scheme. This is a result of the fact that there is more non-zero terms in the global matrix in the former case (at least 1.5 times more non-zero terms). Albeit this increase in the number of non-zero terms in the global matrix, the construction of the L and U factors for the upwind DG discrete operator does not require more CPU than that of the corresponding factors for the centered DG discrete operator. However, a difference in CPU times is clearly observed in the solution time.

Disk-shaped scatterer

We consider now the test problem of the scattering of a plane wave by an infinite elastic cylinder. The analytical solution of this problem is given in annex A. This test problem is a little bit more complicated than the previous one because of its geometry. The computational domain Ω is a ring shape between two circular boundaries: the inner circle with radius a = 2000 m which is associated to the free surface boundary Γ a on which we apply the first condition of (2.2), and the outer circle with radius b = 8000 m corresponding to the boundary Γ b which is assumed to be an artificial boundary and on which we apply the second condition of (2.2). The computational domain Ω is represented on fig. 4.15. The features of the homogeneous material are a mass density ρ = 1.0 10 3 kg.m -3 , Lamé's coefficients λ = 8 M P a and µ = 4 M P a, which imply a velocity v p equal to 4.0 10 3 m.s -1 and a S-waves velocity v s equal to 2.0 10 3 m.s -1 . The distance between circles the two of radius a and b such that it corresponds to 1.5 times the wavelength λ o . We remind that λ o = v p f , where f is the frequency. f is chosen equal to 4 Hz. We discretize the computational domain Ω into three unstructured meshes with respectively 1200, 5100 and 21 000 elements. Two of these meshes are shown on figs. [START_REF] Delcourte | A high order discontinuous Galerkin method for the seismic wave propagation[END_REF].16 and 4.17; their characteristics are given in table 4.6. The numerical convergence is presented on fig. 4.18 for the centered DG and upwind DG formulations. We remark that we do not obtain the expected convergence orders for both schemes, and that the two schemes behave similarly. This is probably due to the fact that the geometric error dominates. Indeed, for this test problem, the curved boundaries are discretized by affine elements which a limitation for obtaining higher convergence orders. 20 and4.21 we plot the numerical solution for V x component obtained respectively with the centered and upwind DGFD-P 2 formulations on mesh M2 (with 5100 elements). We show the exact solution on fig. 4.19. It is interesting to look at the absolute error between the numerical solution and the exact solution (fig. 4.22 for the centered scheme and fig. 4.23 for the upwind scheme). As for the propagation of the plane wave in the homogeneous medium, we plot a 1D cut of the numerical solution (green cross), figs. [START_REF] Delcourte | A high order discontinuous Galerkin method for the seismic wave propagation[END_REF].24 and 4.25 and we compare it to the exact solution (red line). The computational performances of both methods are summarized in tables 4.7 for the number of non-zero terms in the global matrix and the memory used (in MB), 4.8 for the time required for the construction of the global matrix (in seconds) and for the resolution (in seconds) and 4.9-4.10 for the mean and relative errors on V x and σ xx components.

Mesh # Mesh elements # Mesh vertices

Elastic solid scatterer problem

Finally we consider the problem of the scattering of a plane wave by an elastic solid in an infinite circle, which corresponds to an heterogeneous wave propagation problem. The analytical solution of this problem is given in annex B. The computational domain Ω = Ω a ∪Ω b is represented The numerical convergence is presented on fig. 4.29 for the centered DG and upwind DG formulations. We obtain the same numerical convergence for both methods as with the previous test problem.

Inria # Mesh
On figs. [START_REF] Delcourte | A high order discontinuous Galerkin method for the seismic wave propagation[END_REF].31 and 4.32 we plot the numerical solution for V x component obtained respectively with the centered and upwind DGFD-P 2 formulations on mesh M2 (with 5100 elements). We show the exact solution on fig. The computational performances of both methods are summarized in tables 4.12 for the number of non-zero terms in the global matrix and the memory used (in MB), 4.13 for the time required for the construction of the global matrix (in seconds) and for the resolution (in seconds) and 4.14-4.15 for the mean and relative errors on V x and σ xx components.

Inria 5 Conclusion

In this report, we studied two nodal DG methods, the centered flux DG method and the upwind flux DG method, for the solution of the 2D elastic Helmholtz equations. We numerically observed that the centered DG formulation is suboptimal in terms of convergence speed as compared to a classical FE method or the upwind DG formulation considered here. Moreover accuracy results are generally better with the upwind flux DG method for a given mesh. Indeed, the upwind DG method is more accurate than the centered DG one whatever is the computational domain configuration or the medium complexity.

However the upwind DG scheme required more memory space (around 1,5 times much memory) and computational time (around twice much time). This is due to the fact that, in the upwind flux case, the number of non-zero terms in the discretematric operator is around twice as important as the one of the centered method.

The facts that the centered DG method is not enough accurate and the upwind DG method is expensive in terms of memory and computational time, motivate the study of an alternative DG formulation: the hybridizable DG (HDG) method. In a subsequent report, we present and study numerically such a HDG method the 2D elastic Helmholtz equations and we compare the HDG results with those obtained with upwind DG method which we consider as a reference method for our study.

A Analytical solution of the disk-shaped scatterer problem

In this section, we recall the analytical expression of the solution for the problem of the scattering of a plane wave by an elastic disk-shaped configuration. This analytical solution is expressed in the form of Fourier series. In the case of an infinite solid domain, the total displacement field u can be expressed using two others displacement fields u 1 and u 2 as

u = u 1 + u 2 . (A.1)
Each of these displacement fields u j , j = 1, 2, can be written with the help of two potentials φ j and ψ j u 1 = ∇φ 1 + (-e z ) × ∇ψ 1 , (A.2)

u 2 = ∇φ 2 + (-e z ) × ∇ψ 2 , (A.3)
where

φ 1 = +∞ n=0 A 1 n H (1) n (k p r) cos(nθ), ψ 1 = +∞ n=0 A 2 n H (1) n (k s r) sin(nθ), φ 2 = +∞ n=0 A 3 n H (2) n (k p r) cos(nθ), ψ 2 = +∞ n=0 A 4 n H (2) n (k s r) sin(nθ),
and e z is the third vector of the cartesian basis. H (1) n and H (2) n respectively represent Hankel's functions of first and second kind, defined such as H (1) n (x) = J n (x) + iY n (x) H (2) n (x) = J n (x) -iY n (x). This allows us to write u j = u j r e r + u j θ e θ , (A.5) Since the polar basis vectors e r and e θ are given in the cartesian basis by (cos θ, sin θ) t and (-sin θ, cos θ) t respectively, it follows that the components of the displacement field can be expressed in the cartesian basis as

u = u 1 r + u 2 r e r + u 1 θ + u 2 θ e θ , (
u x = u r cos θ -u θ sin θ = (u 1 r + u 2 r ) cos θ -(u 1 θ + u 2 θ ) sin θ, u z = u r sin θ + u θ cos θ = (u 1 r + u 2 r ) sin θ + (u 1 θ + u 2 θ ) cos θ, (A.7) RR n°8989
with

u 1 r = +∞ n=0 A 1 n k p H (1) n (k p r) + n r A 2 n H (1) n (k s r) cos(nθ), u 1 θ = +∞ n=0 -A 1 n n r H (1) n (k p r) -A 2 n k s H (1) n (k s r) sin(nθ), u 2 r = +∞ n=0 A 3 n k p H (2) n (k p r) + n r A 4 n H (2) n (k s r) cos(nθ), u 2 θ = +∞ n=0 -A 3 n n r H (2) n (k p r) -A 4 n k s H (2) n (k s r) sin(nθ).
Knowing that derivatives of Hankel's functions are determined by one of the two recurrence relations

H (j) n (x) = -H (j) n+1 (x) for n = 0, H (j) n-1 (x) - n x H (j)
n (x) for n > 0, or

H (j) n (x) = -H (j) n+1 (x) for n = 0, -H (j) n+1 (x) + n x H (j) n (x) for n > 0,
we choose to write arbitrarily

H (1) n (x) = -H (1) n+1 (x)
for n = 0, -H

n+1 (x) + n x H (1) n (x) for n > 0, (1)
H (2) n (x) = -H and
n+1 (x) for n = 0, H

n-1 (x) -n x H (2) n (x) for n > 0.

Finally this leads to

u r = +∞ n=0 A 1 n k p H (1) n (k p r) + A 2 n n r H (1) n (k s r) + A 3 n k p H (2) n (k p r) + A 4 n n r H (2) n (k s r) cos(nθ) u θ = +∞ n=0 -A 1 n n r H (1) n (k p r) -A 2 n k s H (1) n (k s r) -A 3 n n r H (2) n (k p r) -A 4 n k s H (2) n (k s r) sin(nθ). (A.8)
In order to determine coefficients A 1 n , A 2 n , A 3 n and A 4 n , we have to use the boundary conditions. We get u = -u inc on Γ a , (A.9)

σn = v p (v • n) n + v s (v • t) t on Γ b . (A.10)
We recall the reader to refer to figure 4.15 for the definition of the boundaries Γ a and Γ b . u inc is the incident wave which is written as

u inc = ∇φ inc + (-e z ) × ∇ψ inc ,

Inria

with ψ inc = 0 and

φ inc = +∞ n=0 ε n i n J n (k p r) cos(nθ)
where

ε n = 1 if n = 0, 2 if n ≥ 1,
. And finally we write

u inc = ∂φ inc ∂r e r + 1 r ∂φ inc ∂θ e θ = +∞ n=0 ε n i n cos(nθ)k p J n (k p r)e r + +∞ n=0 -n r ε n i n J n (k p r) sin(nθ)e θ .
(A.11)

If we develop the equation

σn = v p (v • n) n + v s (v • t) t,
in polar coordinates, we find

σ rr = iωρv p u r , σ rθ = iωρv s u θ . (A.12)
σ rr and σ rθ can be expressed in terms of potentials ψ j et φ j

σ rr = σ 1 rr + σ 2 rr = j=1,2 σ φ j rr + σ ψ j rr , σ rθ = σ 1 rθ + σ 2 rθ = j=1,2 σ φ j rθ + σ ψ j rθ , (A.13)
where

σ φ j rr = λ∆φ j + 2µ ∂ 2 φ j ∂r 2 , σ ψ j rr = 2µ ∂ ∂r 1 r ∂ψ j ∂θ , σ φ j rθ = 2µ 1 r ∂ 2 φ j ∂θ∂r - 1 r 2 ∂φ j ∂θ , σ ψ j rθ = µ 1 r 2 ∂ 2 ψ j ∂θ 2 -r ∂ ∂r 1 r ∂ψ j ∂r .
where φ j satisfies the Helmholtz equation

∆φ j = -k 2 p φ j .
This yields

∆φ j = -k 2 p +∞ n=0 A h n H (j) n (k p r) cos(nθ), h = 1 or 3.
Then, by some calculations, we obtain

∂ 2 φ j ∂r 2 = +∞ n=0 A h n k 2 p H (j) n (k p r) cos(nθ), ∂ ∂r 1 r ∂ψ j ∂θ = - 1 r 2 ∂ψ j ∂θ + 1 r ∂ 2 ψ j ∂r∂θ = - 1 r 2 +∞ n=0 A l n H (j) n (k s r)(n cos(nθ)) + 1 r +∞ n=0 nA l n k s H (j) n (k s r) cos(nθ), l = 2 ou 4, ∂ 2 φ j ∂θ∂r = +∞ n=0 A h n k p H (j) n (k p r) (-n sin(nθ)) , ∂φ j ∂θ = +∞ n=0 A h n H (j) n (k p r) (-n sin(nθ)) , ∂ 2 ψ j ∂θ 2 = +∞ n=0 A l n H (j) n (k s r)(-n 2 sin(nθ)), ∂ ∂r 1 r ∂ψ j ∂r = - 1 r 2 ∂ψ j ∂r + 1 r ∂ 2 ψ j ∂r 2 = - 1 r 2 +∞ n=0 A l n k s H (j) n (k s r) sin(nθ) + 1 r +∞ n=0 A l n k 2 s H (j) n (k s r) sin(nθ).
In summary, for σ 1 , we get

                           σ 1 rr = +∞ n=0 A 1 n k 2 p -λH (1) n (k p r) + 2µH (1) n (k p r) + A 2 n 2µn r - 1 r H (1) n (k s r) + k s H (1) n (k s r) cos(nθ), σ 1 rθ = µ +∞ n=0 -2A 1 n n r k p H (1) n (k p r) - 1 r H (1) n (k p r) + A 2 n - n 2 r 2 H (1) n (k s r) + 1 r k s H (1) n (k s r) -k 2 s H (1) n (k s r) sin(nθ). (A.14) In addition, since k 2 p (λ + 2µ) = k 2 s µ, we obtain λ = k 2 s k 2 p µ -2µ
. Moreover, we remark that

H (1) n (k p r) = - 1 k p r H (1) n (k p r) -(1 - n 2 (k p r) 2 )H (1) n (k p r) = 1 k p r H (1) n+1 (k p r) - n (k p r) 2 H (1) n (k p r) -(1 - n 2 (k p r) 2 )H (1) n (k p r). Inria with C n (:, 1) =          k p H (1) n (k p a) - n a H (1) n (k p a) n 2 -n - 1 2 k 2 p b 2 H (1) n (k p b) + k p bH (1) n+1 (k p b) -iω 2 ρ b 2 2µ H (1) n (k p b) -n (n -1) H (1) n (k p b) -k p bH (1) n+1 (k p b) + iωρv s n b 2µ H (1) n (k p b)          , C n (:, 2) =          n a H (1) n (k s a) -k s H (1) n (k s a) n (n -1) H (1) n (k s b) -k s bH (1) n+1 (k p b) -iωρv p n b 2µ H (1) n (k s b) - n 2 -n - 1 2 k 2 s b 2 H (1) n (k s b) + k s bH (1) n+1 (k s b) + iω 2 ρ b 2 2µ H (1) n (k s b)          , C n (:, 3) =          k p H (2) n (k p a) - n a H (2) n (k p a) -k 2 s -2k 2 p H (2) n (k p b) + 2k 2 p H (2) n (k p b) -iω 2 ρH (2) n (k p b) b 2 2 -nb k p H (2) n (k p b) - 1 b H (2) n (k p b) + iωρv s n b 2µ H (2) n (k p b)          , C n (:, 4) =          n a H (2) n (k s a) -k s H (2) n (k s a) nb -1 b H (2) n (k s b) + k s H (2) n (k s b) -iωρv p n b 2µ H (2) n (k s b) - n 2 b 2 H (2) n (k s b) + 1 b k s H (2) n (k s b) -k 2 s H (2) n (k s b) + iω 2 ρ 1 µ H (2) n (k s b) b 2 2          , A n =     A 1 n A 2 n A 3 n A 4 n     , B n =      -ε n i n k p J n (k p a) n a ε n i n J n (k p a) 0 0     
.

B Analytical solution of the elastic solid scatterer problem

Like in the previous section, we recall here the analytical solution in the form of Fourier series of the elastic disk-shaped solid scatterer problem. The displacement field u is represented with the help to two potentials φ and ψ u = ∇φ + (-e z ) × ∇ψ, (B.1)

where, in the case of a circle, like the domain Ω a , we refer to figure 4.26 for the domains Ω, Ω a , Ω b

φ = φ 3 = +∞ n=0 A 3 n J n (k p r) cos(nθ), ψ = ψ 3 = +∞ n=0 A 4 n J n (k s r) sin(nθ),
and in the case of a ring like

Ω b φ = φ 1 + φ 2 = +∞ n=0 A 1 n H (1) n (k p r) cos(nθ) + +∞ n=0 A 5 n H (2) n (k p r) cos(nθ), ψ = ψ 1 + ψ 2 = +∞ n=0 A 2 n H (1) n (k s r) sin(nθ) + +∞ n=0 A 6 n H (2) n (k s r) sin(nθ).
e z is the third vector of the cartesian basis. 1) n represents the Hankel's function of first kind defined such as H (1) n (x) = J n (x) + iY n (x), (B.2)

k p = ω v p et k s = ω v s . H (
with J n et Y n Bessel's functions respectively of first and second kind. H (2) n is the Hankel's function of second kind defined such as H (2) n (x) = J n (x) -iY n (x). Since polar vectors e r and e θ are defined in cartesian basis by (cos θ, sin θ) t and (-sin θ, cos θ) t respectively, we obtain the components of the displacement field in the cartesian basis

u x = u r cos θ -u θ sin θ, u z = u r sin θ + u θ cos θ. (B.5)

Inria

The free surface condition on Γ a is given by Since C n is inversible, we have to solve the following system to find coefficients A j n , j = 1, 6

σ rr12 + σ inc rr = σ rr3 , σ rθ12 + σ inc rθ = σ rθ3 , (
A n = C -1 n B n . (B.27)

  ) with Γ a ∪ Γ b = ∂Ω and Γ a ∩ Γ b = ∅. The first relation of (2.2) describes a free surface condition on Γ a while the second relation indicates an absorbing boundary condition on Γ b that we will detail in section 3.1.4.

Figure 4 . 1 :

 41 Figure 4.1: Discretization of Ω: mesh M1, 3000 elements.

Figure 4 . 2 :Table 4 . 1 :

 4241 Figure 4.2: Discretization of Ω: mesh M2, 10000 elements.

Figure 4 . 3 :

 43 Figure 4.3: Exact solution, V x component.

Figure 4 . 4 :

 44 Figure 4.4: Numerical solution, mesh M1, centered flux DGFD-P 1 cheme, V x component.

Figure 4 . 5 :

 45 Figure 4.5: Numerical solution, mesh M1, upwind flux DGFD-P 1 scheme, V x component.

Figure 4 . 6 :

 46 Figure 4.6: Numerical solution, mesh M1, centered flux DGFD-P 2 scheme, V x component.

Figure 4 . 7 :

 47 Figure 4.7: Numerical solution, mesh M1, upwind flux DGFD-P 2 scheme, V x component.

Figure 4 . 8 :

 48 Figure 4.8: Numerical solution, mesh M2, centered flux DGFD-P 1 scheme, V x component.

Figure 4 .

 4 Figure 4.10: x-wise distribution of V x , mesh M1: in red line the exact solution, in green cross the centered (left) and upwind (right) DGFD-P 1 solution.

Figure 4 .

 4 Figure 4.11: x-wise distribution of V x , mesh M1: in red line the exact solution, in green cross the centered (left) and upwind (right) DGFD-P 2 solution.

Figure 4 .

 4 Figure 4.12: x-wise distribution of V x , mesh M2: in red line the exact solution, in green cross the centered (left) and upwind (right) DGFD-P 1 solution.

Figure 4 .

 4 Figure 4.13: x-wise distribution of V x , mesh M2: in red line the exact solution, in green cross the centered (left) and upwind (right) DGFD-P 1 solution.

Figure 4 .

 4 Figure 4.14: Convergence order of the centered (left) and upwind (right) flux DGFD methods for plane wave propagation in a homogeneous medium.

Figure 4 .

 4 Figure 4.15: Configuration of the computational domain Ω for the elastic disk-shaped scatterer.

Figure 4 .

 4 Figure 4.16: Discretization of the computational domain Ω for the elastic disk-shaped scatterer: mesh M1, 1200 elements.

Figure 4 .

 4 Figure 4.17: Discretization of the computational domain Ω for the elastic disk-shaped scatterer: mesh M2, 5100 elements.

Figure 4 .

 4 Figure 4.18: Convergence order of the centered (left) and upwind (right) flux DGFD methods for the elastic disk-shaped scatterer.

Figure 4 .

 4 Figure 4.19: Exact solution of V x for the elastic disk-shaped scatterer.

Figure 4 .

 4 Figure 4.20: Numerical solution of V x with the centered DGFD-P 2 method for the elastic disk-shaped scatterer.

Figure 4 .

 4 Figure 4.21: Numerical solution of V x with the upwindd DGFD-P 2 method for the elastic disk-shaped scatterer.

Figure 4 .

 4 Figure 4.22: Absolute error between the exact solution and the solution computed with the centered DGFD-P 2 method for the elastic disk-shaped scatterer.

Figure 4 .

 4 Figure 4.23: Absolute error between the exact solution and the solution computed with the upwind DGFD-P 2 metho for the elastic disk-shaped scatterer.

Figure 4 .

 4 Figure 4.24: x-wise distribution V x : centered DGFD-P 2 solution for the elastic disk-shaped scatterer.

Figure 4 .

 4 Figure 4.25: x-wise distribution of V x : upwind DGFD-P 2 solution for the elastic disk-shaped scatterer.

Inria on fig. 4 .Figure 4 .

 44 Figure 4.26: Configuration of the computational domain Ω for the elastic solid scatterer.

Figure 4 .

 4 Figure 4.27: Discretization of the computational domain Ω for the elastic solid scatterer: mesh M1, 1300 elements.

Figure 4 .

 4 Figure 4.28: Discretization of the computational domain Ω for the elastic solid scatterer: mesh M2, 5400 elements.

Figure 4 .

 4 Figure 4.29: Convergence order of the centered (left) and upwind (right) flux DGFD methods for the elastic solid scatterer.

Figure 4 .

 4 Figure 4.31: Numerical solution of V x with the centered DGFD-P 2 method for the elastic solid scatterer.

Figure 4 .Figure 4 .

 44 Figure 4.32: Numerical solution of V x with the upwind DGFD-P 2 method for the elastic solid scatterer.

Figure 4 .

 4 Figure 4.33: Absolute error between the exact solution and the solution computed with the centered DGFD-P 2 method for the elastic solid scatterer.

Figure 4 .

 4 Figure 4.34: Absolute error between the exact solution and the solution computed with the upwind DGFD-P 2 method for the elastic solid scatterer.

Figure 4 .

 4 Figure 4.35: x-wise distribution V x : centered DGFD-P 2 solution for the elastic solid scatterer.

Figure 4 .

 4 Figure 4.36: x-wise distribution of V x : upwind DGFD-P 2 solution for the elastic solid scatterer.

(A. 4 )

 4 The Hankel function of second kind is the conjugate of the Hankel function of first kind. J n and Y n are respectively Bessel's functions of first and second kind. k p = ω v p is the P -wave number and k s = ω v s the S-wave number. In polar coordinates, we have ∇φ j = ∂φ j ∂r e

  u = u r e r + u θ e θ ,

n 1 rn

 1 B.15) and the absorbing condition on Γ b byσ rr12 = iωv p u 12r , σ rθ12 = iωv s u 12θ . (B.16)The components σ rr and σ rθ of σ are expressed such asσ rr = σ φ rr + σ ψ rr , σ rθ = σ φ rθ + σ ψ rθ , C n (k p r) cos(nθ), h = 1, 3 or5, and C n (k p r) = H(1)n (k p r), J n (k p r) or H(2)n (k p r). For σ 1 we get  r) 2 n + (k p r) 2 -n 2 H (1) n (k p r) (k s r) -k s H (1) n+1 (k s r) cos(nθ), σ rθ1 = µ -n 2 + n + 1 2 (k s r) 2 H (1) n (k s r) -k s H (1) n+1 (k s r) sin(nθ). r) 2 n -(k p r) 2 + n 2 H (2) n (k p r) -(k s r) + k s H (2)n-1 (k s r) cos(nθ).

r n 2 + n -1 2 (k s r) 2 H ( 2 ) 1 r ( 1 - 4 n 2 r 1 rε

 22221141 n (k s r) + k s H (2) n-1 (k s r) sin(nθ). (B.19) r) 2 n + (k p r) 2 -n 2 J n (k p r) + 2µ 1 k p r J n+1 (k p r) n) J n (k s r) -k s J n+1 (k s r) cos(nθ), J n+1 (k p r) + n -1 r J n (k p r) + A -n 2 + n + 1 2 (k s r) 2 J n (k s r) -k s J n+1 (k s r) sin(nθ). n i n -λk 2 p -2µ 1 r 2 nn 2 + (k p r) 2 J n (k p r) + 2µ k p r J n+1 (k p r) cos(nθ), σ inc rθ = +∞ n=0 -2µε n i n n r -1 -n r J n (k p r) -k p J n+1 (k p r) sin(nθ).(B.21) For r = a we have+∞ n=0 A 1 n k p -H

k 2 s µ, we replace λ by k 2 s k 2 pµε

 2 -2µ and rewriting equations (B.22) to (B.25)under matrix form, we finally getC n A n = B n , n i n k p J n (k p a) -n a ε n i n J n (k p a) ε n i n k 2 s J n (k p r) n (k p r) -n 2 r 2 J n (k p r) -2ε n i n n r k p J n (k p r) -1 r J n (k p r)

Table 4 .

 4 2: Number of non-zero terms in the global matrix and memory used.

	# Mesh elements Interpolation	Non-zeros terms	Memory (MB)
		degree
			Cent. scheme Upw. scheme Cent. scheme Upw. scheme
	3100	1	7.5e+05	1.5e+06	204	288
	10300	1	2.5e+06	5.1e+06	877	1076
	45000	1	1.1e+07	2.2e+07	4489	5492
	3100	2	2.2e+06	4.3e+06	527	804
	10300	2	7.4e+06	1.4e+07	2036	3097
	45000	2	3.2e+07	6.2e+07	10687	15965
	3100	3	5.8e+06	9.4e+06	1246	1656
	10300	3	1.9e+07	3.1e+07	5020	6600
	45000	3	8.3e+07	1.3e+08	27228	34597
	3100	4	1.2e+07	1.8e+07	1980	2749
	10300	4	4.1e+07	5.9e+07	7372	10098
	45000	4	1.8e+08	2.6e+08	37018	50297
						Inria

Table 4 .

 4

	h (m) Interpolation	Mean Error V x	Relative Error V x
	degree	Cent. scheme Upw. scheme Cent. scheme Upw. scheme
	625.0	4.8e-02	1.8e-02	15.9	6.0
	312.5	3.2e-03	1.7e-03	3.3	1.8
	56.25	1.7e-04	9.9e-05	0.8	0.4
	625.0	1.7e-03	1.6e-03	0.5	0.5
	312.5	5.9e-05	5.7e-05	6.0e-02	5.8e-02
	56.25	1.4e-06	1.4e-06	6.1e-03	6.1e-03
	625.0	1.4e-04	1.3e-04	4.2e-02	4.1e-02
	312.5	3.2e-06	3.1e-06	3.3e-03	3.2e-03
	56.25	3.8e-08	3.7e-08	1.7e-04	1.7e-04
	625.0	1.1e-05	1.0e-05	3.4e-03	3.2e-03
	312.5	7.6e-08	6.8e-08	7.8e-05	7.0e-05
	56.25	3.5e-10	3.1e-10	1.6e-06	1.4e-06

3: Time required for the global matrix construction and for the system resolution.

Table 4 .

 4 4: Mean and relative errors on V x .

	RR n°8989

Table 4 .

 4

5: Mean and relative errors on σ xx and convergence order.

Table 4 . 6

 46

	h min	h max	h min /h max

: Characteristics of the three meshes On figs. 4.

Table 4 .

 4 7: Number of non-zero terms in the global matrix and memory used.

	elements Interpolation	Non-zeros terms	Memory (MB)
		degree
			Cent. scheme Upw. scheme Cent. scheme Upw. scheme
	1200	1	2.9e+05	5.9e+05	72	104
	5100	1	1.2e+06	2.5e+06	379	496
	21000	1	4.9e+06	1.0e+07	1895	2292
	1200	2	8.6e+05	1.6e+06	180	269
	5100	2	3.7e+06	7.0e+06	925	1360
	21000	2	1.5e+07	2.7e+07	4361	6578
	1200	3	2.2e+06	3.6e+06	391	525
	5100	3	9.4e+06	1.5e+07	2154	2921
	21000	3	3.8e+07	6.2e+07	11033	14131
	1200	4	4.7e+06	6.8e+06	673	895
	5100	4	2.0e+07	2.9e+07	3393	4537
	21000	4	8.1e+07	1.2e+08	15679	21186
	# Mesh elements Interpolation	Construction time (s)	Solution time (s)
		degree
			Cent. scheme Upw. scheme Cent. scheme Upw. scheme
	1200	1	1.1e-02	2.5e-02	0.4	0.7
	5100	1	4.8e-02	0.1	2.5	4.0
	21000	1	0.2	0.4	16.0	25.1
	1200	2	3.5e-02	7.1e-2	1.3	2.6
	5100	2	0.1	0.3	7.8	14.7
	21000	2	0.6	1.2	47.7	93.7
	1200	3	8.8e-02	0.2	2.9	5.3
	5100	3	0.4	0.7	18.2	38.2
	21000	3	1.5	2.7	129.4	249.3
	1200	4	0.2	0.3	5.5	10.2
	5100	4	0.8	1.4	33.6	65.6
	21000	4	3.3	5.5	234.5	447.3

Table 4 .

 4 8: Time required for the global matrix construction and for the system resolution.

	RR n°8989

Table 4 .

 4 9: Mean and relative errors on V x .

	h (m) Interpolation Mean Error σ xx Relative Error σ xx Convergence order
	degree	Cent.	Upw.	Cent.	Upw.	Cent.	Upw.
		scheme scheme scheme	scheme	scheme	scheme
	1016.9	5950.9	3999.1	294.6	198.0	-	-
	490.1	1001.7	829.2	166.1	137.5	0.3	0.2
	245.9	133.3	66.1	74.4	36.9	0.9	1.5
	1016.9	5673.7	3802.5	193.9	128.0	-	-
	490.1	448.3	159.8	58.8	20.9	1.3	2.0
	245.9	6.5	4.8	3.4	2.5	4.3	3.2
	1016.9	2372.0	1195.1	72.0	36.3	-	-
	490.1	56.5	37.9	7.3	4.9	3.1	2.7
	245.9	3.1	3.0	1.6	1.5	2.5	1.9
	1016.9	1141.6	743.8	34.3	22.3	-	-
	490.1	31.6	31.0	4.1	3.4	3.0	2.4
	245.9	3.0	3.0	1.5	1.5	1.6	1.6

Table 4 .

 4 10: Mean and relative errors on σ xx and convergence order.

Table 4 .

 4 4.30. We plot the absolute error between the numerical solution Mesh # Mesh elements # Mesh vertices h min h max h min /h max 11: Characteristics of the three meshes

	M1	1300	700	440.5 1016.9	2.3
	M2	5400	2800	211.6	490.1	2.3
	M3	22 000	11 000	105.9	245.9	2.3

Table 4 .

 4 12: Number of non-zero terms in the global matrix and memory used.

	Inria

Table 4 .

 4 13: Time required for the global matrix construction and for the system resolution.

	RR n°8989

Table 4 .

 4 14: Mean and relative errors on V x .

	h (m) Interpolation Mean Error σ xx Relative Error σ xx Convergence order
	degree	Cent.	Upw.	Cent.	Upw.	Cent.	Upw.
		scheme scheme scheme	scheme	scheme	scheme
	1016.9	6482.2	4119.1	261.3	166.1	-	-
	490.1	1125.5	830.8	154.0	113.7	0.3	0.2
	245.9	125.8	63.9	60.2	30.6	1.1	1.6
	1016.9	5934.8	3742.9	172.1	108.6	-	-
	490.1	459.3	143.2	52.5	16.4	1.4	2.1
	245.9	6.3	5.5	2.9	2.5	4.3	2.9
	1016.9	2353.6	1054.3	62.4	27.9	-	-
	490.1	56.0	41.6	6.3	4.7	3.2	2.4
	245.9	3.9	3.9	1.7	1.7	2.1	1.6
	1016.9	1073.5	667.3	28.3	17.6	-	-
	490.1	26.4	26.9	3.0	3.0	3.3	2.6
	245.9	3.8	3.8	1.7	1.7	0.9	0.9

Table 4 .

 4 15: Mean and relative errors on σ xx and convergence order.

  J n (k p r) sin(nθ) (k s r) + k s J n (k s r) -(k s r) + k s H (1) n (k s r) cos(nθ) = ε n i n (λ + 2µ) k 2 p J n (k p r) + λ 1 r k p J n (k p r) -n 2 r 2 J n (k p r) cos(nθ). J n (k s r) -k 2 s J n (k s r) +

		A 3 n (λ + 2µ) k 2 p J n (k p r) + λ	1 r	k p J n (k p r) -	n 2 r 2 J n (k p r)	+
		2µA 4 n J n A 1 n r -1 r n (λ + 2µ) k 2 p H (1) n	(k p r) + λ	1 r	k p H (1) n (k p r) -	n 2 r 2 H (1) n (k p r)	-	(B.24)
		2µA 2 n	n r	-1 r	H (1)
				-2µA 3 n	n r	k p J n (k p r) -	1 r	J n (k p r) +
					µA 4 n -k s 2µA 1 n 2 r 2 J n (k s r) + 1 r n n r k p H (1) n (k p r) -1 r	H (1) n (k p r) -
					µA 2 n -	n 2 r 2 H (1) n (k s r) +	1 r	k s H (1)
								n+1 (k p r) +	n k p r	H (1) n (k p r) +	n r	A 2 n H (1) n (k s r)+
						A 5 n k p H n-1 (k p r) -(2)	n k p r	H (2) n (k p r) +	n r	A 6 n H (2) n (k s r)-
	A 3 n k p -J n+1 (k p r) + n J +∞ n k p r J n (k p r) -n r A 4 n=0 -A 1 n n r H (1) n (k p r) -A 2 n k s -H (1) n+1 (k s r) + n k s r H (1) n (k s r) -A 5 n	n r	H (2) n (k p r)-
			A 6 n k s H	(2) n-1 (k s r) -	n k s r	H (2) n (k s r) + A 3 n	n r	J n (k p r)+
			A 4 n k s -J n+1 (k s r) +	n k s r	J n (k s r) sin(nθ)
	=	+∞ n=0	n r	ε n i n Inria

n (k s r) cos(nθ),

= +∞ n=0 -ε n i n k p -J n+1 (k p r) + n k p r J n (k p r) cos(nθ), n n (k s r) -k 2 s H (1) n (k s r) sin(nθ) = -2µε n i n n r k p J n (k p r) -1 r J n (k p r) sin(nθ).

RR n°8989

avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements

The authors acknowledge the support by the INRIA-TOTAL strategic action DIP (dip.inria.fr).

It follows that

n+1 (k p r) cos(nθ),

n+1 (k p r) -

n+1 (k s r) sin(nθ).

(A.15)

We can do the same for σ 2 and we obtain

If we sum these two systems, we finally get for σ

(A.17)

To summarize, using boundary conditions, we get

Since n represents modes of Fourier's serie, we compute the coefficients A j n , j = 1, 4 by solving the following system at each Fourier mode n, since C n is inversible

where C n = C n (:, 1) C n (:, 2) C n (:, 3) C n (:, 4)

Inria

Then, in Ω a , the components of displacement field u 3 are written

where

In Ω b , the components of displacement field u 1 are written

where

so we get

Derivatives of Hankel's function are determined by the following reccurence relations

n+1 (x) + n x H (1) n (x) for n > 0, (2) n (x) for n > 0.

In the same way, the derivative of the Bessel's function of first kind is

Using these definitions we rewrite u ir and u iθ , i = 12, 3

n-1 (k p r) -

In order to determine coefficients A i n , i = 1, .., 6, we use boundary conditions. We get

on Γ a , (B.12)

We refer to figure 4.26 for the definition of boundaries Γ a and Γ b .

u inc is the displacement field of the incident wave described by

with ψ inc = 0 and

where

In summary, we have

n-1 (k s r) -

n-1 (k p r) -

n-1 (k s r) + 2A


Boyan Kolev

Raquel Pau

Oleksandra Levchenko

Patrick Valduriez

Ricardo Jiménez-Peris

José Pereira

Benchmarking Polystores: the CloudMdsQL Experience

Keywords: benchmark, cloud, heterogeneous data, polystore

The CloudMdsQL polystore provides integrated access to multiple heterogeneous data stores, such as RDBMS, NoSQL or even HDFS through a big data analytics framework such as MapReduce or Spark. The CloudMdsQL language is a functional SQL-like query language with a flexible nested data model. A major capability is to exploit the full power of each of the underlying data stores by allowing native queries to be expressed as functions and involved in SQL statements. The CloudMdsQL polystore has been validated with a good number of different data stores: HDFS, key-value, document, graph, RDBMS and OLAP engine. In this paper, we introduce the benchmarking of the CloudMdsQL polystore and evaluate the performance benefits of important features enabled by the query language and engine.

I. INTRODUCTION

The blooming of different cloud data management infrastructures, specialized for different kinds of data and tasks, has led to a wide diversification of DBMS interfaces and the loss of a common programming paradigm. This makes it very hard for a user to integrate and analyze her data sitting in different data stores, e.g. RDBMS, NoSQL, and HDFS. The CoherentPaaS project [1] addresses this problem, by providing a rich platform integrating different data management systems specialized for particular tasks, data and workloads. The platform is designed to provide a common programming model and language to query multiple data stores.

The problem of accessing heterogeneous data sources has long been studied in the context of multidatabase and data integration systems [START_REF] Özsu | Principles of Distributed Database Systems[END_REF]. More recently, with the advent of cloud databases and big data processing frameworks, the solution has evolved towards polystores (also called multistore systems) that provide integrated access to a number of RDBMS, NoSQL and HDFS data stores through a common query engine. Data mediation SQL engines, such as Apache Drill, Spark SQL, and SQL++ provide common interfaces that allow different data sources to be plugged in (through the use of wrappers) and queried using SQL. The BigDAWG polystore [START_REF] Duggan | The BigDAWG polystore system[END_REF] goes farther by enabling queries across "islands of information", where each island corresponds to a specific data model and its language and provides transparent access to a subset of the underlying data stores through the island's data model. Another family of multistore systems [START_REF] Dewitt | Split query processing in Polybase[END_REF][START_REF] Lefevre | MISO: souping up big data query processing with a multistore system[END_REF] has been introduced with the goal of tightly integrating big data analytics frameworks (e.g. MapReduce or Spark) with traditional RDBMS, by sacrificing the extensibility with other data sources. However, since none of these approaches supports the ad-hoc usage of native queries, they do not preserve the full expressivity of an arbitrary data store's query language. But what we want to give the user is the ability to express powerful ad-hoc queries that exploit the full power of the different data store languages, e.g. directly express a path traversal in a graph database.

The Cloud Multidatastore Query Language (CloudMdsQL) is a functional SQL-like query language, designed to serve the querying capabilities of the CoherentPaaS platform. A CloudMdsQL query can address multiple heterogeneous databases by means of nested subqueries [START_REF] Kolev | CloudMdsQL: querying heterogeneous cloud data stores with a common language[END_REF]. Each subquery addresses directly a particular data store and may contain embedded invocations to the data store's native query interface. Thus, the major innovation is that a CloudMdsQL query can exploit the full power of local data stores, by simply allowing some local data store native queries (e.g. a breadth-first search query against a graph database) to be called as functions, and at the same time be optimized, e.g. by pushing down select predicates, using bind join, etc.

One of the major challenges in front of the CloudMdsQL language/engine is to allow joins across heterogeneous data stores and to be able to perform them in an efficient way. For this reason, we pay special attention to the use of bind joins [START_REF] Haas | Optimizing queries across diverse data sources[END_REF] and we apply this technique even when native queries are used. In this paper, we introduce the benchmarking of the CloudMdsQL system, which evaluates the ability to run optimized queries across heterogeneous data stores, as well as the functional and performance benefits of the usage of native queries in combination with SQL statements.

The rest of the paper is organized as follows. Section 2 gives an overview of the CloudMdsQL language and engine. Section 3 presents the benchmark environment and the results of the performance evaluation. Section 4 concludes.

II. CLOUDMDSQL OVERVIEW

The CloudMdsQL language is SQL-based with the extended capabilities for embedding subqueries expressed in terms of each data store's native query interface. The common data model respectively is table-based, with support of rich datatypes that can capture a wide range of the underlying data stores' datatypes, such as arrays and JSON objects, in order to handle non-flat and nested data, with basic operators over such composite datatypes.

The query engine follows a mediator/wrapper architecture, where data stores are accessed through wrappers, which implement a common interface. The query compiler decomposes the query into a query execution plan, which appears as a directed acyclic graph of relational operators where leaf nodes correspond to subqueries for the wrappers to execute directly against the data stores.

A. Query Language

Queries that integrate data from several data stores usually consist of subqueries and an integration SELECT statement. A subquery is defined as a named table expression, i.e. an expression that returns a table and has a name and signature. The signature defines the names and types of the columns of the returned relation. Thus, each query, although agnostic to the underlying data stores' schemas, is executed in the context of an ad-hoc schema, formed by all named table expressions within the query. A named table expression can be defined by means of either an SQL SELECT statement (that the query compiler is able to analyze and possibly rewrite) or a native expression (that the query engine considers as a black box and delegates its processing directly to the data store). For example, the following simple CloudMdsQL query contains two subqueries, defined by the named table expressions T1 and T2, and addressed respectively against the data stores rdb (an SQL database) and mongo (a MongoDB database): T1(x int, y int)@rdb = ( SELECT x, y FROM A ) T2(x int, z array)@mongo = {* db.B.find( {$lt: {x, 10}}, {x:1, z:1, _id:0} ) *} SELECT T1.x, T2.z FROM T1, T2 WHERE T1.x = T2.x AND T1.y <= 3

The purpose of this query is to perform relational algebra operations (expressed in the main SELECT statement) on two datasets retrieved from a relational and a document database. The two subqueries are sent independently for execution against their data stores in order the retrieved relations to be joined by the CloudMdsQL query engine. The SQL table expression T1 is defined by an SQL subquery, while T2 is a native expression (identified by the special bracket symbols {* *}) expressed as a native MongoDB call. Note that subqueries to some NoSQL data stores can also be expressed as SQL statements; in such cases, the wrapper must provide the mapping from relational operators to native calls. Within our evaluation, unlike in the example above, we use an SQL wrapper to query MongoDB, which also benefits from subquery rewriting. Note that parametrization and nesting is also available in SQL and native named tables. Within our evaluation, we give an example that involves the Sparksee graph database and we use its Python API to express subqueries that benefit from all of the features described above. In fact, our initial query engine implementation enables Python integration; however support for other languages (e.g. JavaScript) for user-defined operations can be easily added.

B. Bind Join

CloudMdsQL uses bind join as an efficient method for performing semi-joins across heterogeneous data stores that uses subquery rewriting to push the join conditions. For example, the list of distinct values of the join attribute(s), retrieved from the left-hand side subquery, is passed as a filter to the right-hand side subquery. To illustrate it, let us consider the following CloudMdsQL query:

A(id int, x int)@DB1 = (SELECT a.id, a.x FROM a) B(id int, y int)@DB2 = (SELECT b.id, b.y FROM b) SELECT a.x, b.y FROM b JOIN a ON b.id = a.id
Let us assume that the optimizer has decided to use the bind join method and that the join condition will be bound to the right-hand side of the equi-join operation. First, the relation B is retrieved from the corresponding data store using its query mechanism. Then, the distinct values of B.id are used as a filter condition in the query that retrieves the relation A from its data store. Assuming that the distinct values of B.id are b 1 … b n , the query to retrieve the right-hand side relation of the bind join uses the following SQL approach (or its equivalent according to the data store's query language), thus retrieving from A only the rows that match the join criteria: SELECT a.id, a.x FROM a WHERE a.id IN (b1, …, bn)

In the above example, using bind join will be reasonable only in the presence of an index on the column a.id in data store DB1, because such a presence will make the pushed down IN operator avoid an expensive table scan.

The way to do the bind join analogue for native/Python queries is through the use of a JOINED ON clause in the named table signature. For example, if A is defined as the Python function below, as A.id participates in an equi-join, the values b 1 … b n will be provided to the Python code through the iterator b_keys (in this context, we refer to the table B as the "outer" table, and b_keys as the outer keys):

A(id int, x int JOINED ON id REFERENCING OUTER AS b_keys)@DB1 = {*
for id in CloudMdsQL.b_keys: yield ( id, db.get_x(id) ) *}

In fact, each wrapper may provide different mechanisms for the programmer to consume the values of the outer keys. In this example, as the native query is expressed in an imperative language, the set of outer keys is represented by an iterator object; however, for declarative native queries, more appropriate will be to define a placeholder, which the wrapper will substitute with the textual representation of the set of values.

C. Query Engine

For the current implementation of the query engine, we modified the open source Apache Derby database to accept CloudMdsQL queries and transform the corresponding execution plan into Derby SQL operations. We developed the query planner and the query execution controller and linked them to the Derby core, which we use as the operator engine. Derby allows extending the set of SQL operations by means of CREATE FUNCTION statements. This type of statements creates an alias, with an optional set of parameters, to invoke a specific Java component as part of an execution plan. Thus, for each named table expression in a query, a table function is created dynamically, which invokes the corresponding wrapper as a Java class. Thus, Derby handles global execution, delegating local optimization and execution to the underlying data stores. As a second step, the query engine evaluates which named expressions are queried more than once and must be cached into the temporary table storage, which will be always queried and updated from the specified Java functions to reduce the query execution time. Finally, the last step consists of translating all operation nodes that appear in the execution plan into a Derby specific SQL execution plan.

III. BECNHMARKING AND PERFORMANCE EVALUATION

The goal is to assess the ability of the query engine to perform optimized CloudMdsQL queries across data stores. To evaluate the performance improvements, we compare the execution times (with and without optimization) of queries that involve a quite diverse set of data stores. We divide the performed test cases in two groups that focus on the performance benefits thanks to: (1) bind joins and (2) the support of native queries.

A. Data Stores

We performed the experiments on 5 different data stores, each having different interfaces and data models, in order to have a representative data store in each of the categories: relational, key-value, document, graph databases, and an OLAP engine. They are summarized in this subsection, detailing the aspects that are important for the evaluation, namely: (1) the data model, (2) the query language, (3) whether SQL or native queries are used for subquerying, and (4) how the wrapper handles the set of outer keys to provide bind join support.

MonetDB is a column-store database with an SQL engine atop. It follows a relational data model and is queried by SQL through CloudMdsQL subqueries. To handle the outer keys for bind join, the wrapper creates a temporary table (say outer_keys) in the MonetDB database and populates it with the values of the outer keys; then it rewrites the subquery by adding the predicate IN (SELECT * FROM outer_keys). If the bind join is on multiple keys, the generated predicate is different and uses the EXISTS operator.

MongoDB is a document database that has a Java based query interface, where the chain of operations is specified through the use of API calls with JSON documents as parameters that typically return document collections. Since in some particular but very common cases (e.g. collections of flat documents with fixed field names), the document data model can be considered as a subset of the relational model, it is possible to map simple SQL commands to native MongoDB queries. For this reason, our MongoDB wrapper provides a translation from relational operations to native API calls, so that MongoDB can be subqueried by SQL statements. In addition, it also supports native MongoDB queries. When consuming the outer keys for bind join, the wrapper dynamically builds a filter condition using the MongoDB $in operator and applies it through a call to the find() function.

Apache HBase is a key-value data store with a very simple interface mostly consisting of primitives to scan a whole table or to get the value for a particular key. Within CloudMdsQL, subqueries to HBase are expressed with native statements that can be of the form of either (1) scan '<table_name>' or (2) get '<table_name>', <key>. Therefore, there is no way to directly apply an equivalent of an IN operator to HBase subqueries, but the optimized functionality of a bind join can be achieved through the use of scalar lookups.

Sparksee is a graph database that has a Python API and a declarative notation for graph algebra expressions. Within CloudMdsQL, subqueries to Sparksee are expressed in Python code. When handling the outer keys for bind join, the wrapper provides to the Python code an iterator object that the programmer can use to get the values of the outer keys and use them in native API calls to apply the filter.

ActivePivot is an OLAP engine where data are organized in cubes and queried through MDX statements. Subqueries within CloudMdsQL are expressed as native MDX commands. To handle the outer keys for bind join, the programmer uses a special placeholder in the MDX expression, which the wrapper identifies and dynamically replaces with the textual representation of the set of values of the outer keys.

B. Experimental Setup

We performed our experiments using data based on the TPC-H benchmark schema (www.tpc.org/tpch). Each data store contains the same 8 datasets, generated by the TPC-H generator, stored in tables or in the corresponding for the data store structures, as follows. In MonetDB the TPC-H datasets are stored as relational tables, in MongoDB -as document collections. In HBase they are organized in HBase tables, where each key-value pair corresponds to a row from the dataset; in particular, the key corresponds to the value(s) of the key column(s), while the value contains the serialization of the entire row. In Sparksee, each of the TPC-H datasets is loaded into a set of graph nodes of the same type, while relationships between them are mapped to graph edges, i.e., a join between two datasets is expressed as a graph neighborship query. In ActivePivot, all the generated TPC-H data are loaded in one single cube.

For the experiments, we use a cluster of 6 identical machines (8GB RAM, 4 CPU cores @2.4GHz), of which one node is dedicated to the CloudMdsQL engine and one node per each of the data stores. Each data store node contains data generated as per the TPC-H schema, with scale factor 1. Our queries use the datasets LINEITEM and ORDERS, where the columns L_ORDERKEY, L_PARTKEY, O_ORDERKEY, O_ORDERSTATUS and O_ORDERPRIORITY are indexed.

Our queries use the datasets LINEITEM and ORDERS, where the columns L_ORDERKEY, L_PARTKEY, O_ORDERKEY, O_ORDERSTATUS and O_ORDERPRIORITY are indexed.

C. Evaluation of Bind Joins

With this set of test cases, we focus on the evaluation of performance benefits thanks to the use of bind joins, which is an efficient technique to apply semi-joins between datasets from different data stores by allowing the output of one of the subqueries to be used by the other subquery to filter out the unnecessary for the join rows. As we apply this technique to a quite diverse set of data stores, even if they do not natively support SQL querying, we consider this evaluation very important for the benchmarking of the CloudMdsQL system.

We consider the following query, assuming that the column O_ORDERKEY is indexed. The selection predicate L_PARTKEY = 10 results in retrieving only 24 rows from the LINEITEM table which are then supposed to be joined with the ORDERS table, which contains 1.5 million rows. The non-optimized execution of this query implies that the entire ORDERS table is retrieved at the common query engine before being joined with the small table. In the optimized variant, the T2 subquery is rewritten by adding the predicate O_ORDERKEY IN (k 1 , k 2 , …, k 24 ), where k i are the values of the column L_ORDERKEY, taken from the small table T1. Thus, only the rows that match the join condition (24 instead of 1.5 million) will be retrieved from the ORDERS table. The above query, as specified, is applicable only if both datastore1 and datastore2 support SQL subquerying (MonetDB and MongoDB). To apply the bind join method, the query compiler takes care of doing subquery rewriting by pushing down to the subquery for T2 a predicate condition containing an IN operator that takes the values from an intermediate named table, containing the set of outer keys. As described in the previous subsection, the wrappers for MonetDB and MongoDB use different approaches to achieve the goal.

Next, we will explore the performance of this approach when using native subqueries, especially when T2 is native, as it requires the native query to access intermediate data from the table storage of the common query engine. This is done with the help of the programmer by adding to the signature of T2 the following clause:

JOINED ON O_ORDERKEY REFERENCING OUTER AS t1_keys

Thus, whenever T2 is used for a bind join, the join key values of the other side of the join (the outer keys, taken from T1) are provided in the intermediate table t1_keys and the native query for T2 needs to use the corresponding mechanism that its wrapper provides to access these join keys and use them to return only rows that match the join criteria. Sparksee and ActivePivot wrappers provide different mechanisms for expressing this, as illustrated in the corresponding implementations of T2 below (for simplicity, the full signature of T2 is not repeated and some notations are shortened): T2( ... JOINED ON O_ORDERKEY REFERENCING OUTER AS t1_keys)@sparksee = {* rs = CloudMdsQL.t1_keys() while rs.next(): rs2 = graph.compute( \ "GRAPH::SELECT( 'ORDERS'.'O_ORDERKEY' = " + \ rs [1] + ")") while rs2.next(): yield( rs2 [1], rs2 [START_REF] Dewitt | Split query processing in Polybase[END_REF], ..., rs2[9] ) rs2.close(); rs.close() *}

The above Python code for Sparksee iterates through the outer keys and for each key finds the corresponding graph node and yields a tuple with its attributes, which correspond to the columns of the ORDERS dataset. The above approaches allow for the native query to filter out rows by applying an analogue of the IN operator. This is possible, because Sparksee and ActivePivot provide powerful query mechanisms. However, for simpler query notations, such as the one for the key-value data store HBase, this approach is not applicable, as HBase has only primitives either to retrieve a whole table (scan) or to make a lookup for the value for a particular key (get) or ranges of keys. That is why, to achieve the same benefits of bind join as above, when T2 is an HBase query, we take advantage of scalar lookups, that allow a parameterized named table (T2) to be used as a scalar function and evaluated for every value of a column from another table (T1). So, T2 is defined as a parameterized function that gets all the values of ORDERS for a particular key. Then, T2 is called in the SELECT list of the main SELECT statement of the query, instead of being joined with T1. The implementation of T2 for HBase and the main SELECT statement for this case are as follows:

., T2(L_ORDERKEY).O_COMMENT FROM T1

We performed our experiments in the context of the above query and its variations with respect to the T2 subquery for the 5 different data stores. Since the retrieval of the T1 table is quite trivial and its retrieval time is similar for all data stores, we do not focus on various data stores for the T1 subquery and performed all the experiments with MongoDB as datastore1, varying datastore2 among all the 5 data stores.

We also experimented with 5 different selectivity factors of the bind join condition by changing the T1 filter predicate (initially L_PARTKEY = 10), thus varying the number of rows returned by T1, which results in different cardinalities of the set of outer keys. The different filter predicates we used and their corresponding selectivity factors (SF) are detailed in Table 1. Thus, the results of the current or similar benchmark can serve as a statistical base to determine the profitability of bind join compared to traditional join for each data store as a function of the join selectivity, as well as to determine a selectivity threshold, below which using bind join will be considered reasonable. Since each data store and wrapper uses different mechanisms to handle the outer keys, these metrics are quite different for each of the 5 data stores.

Table 2 summarizes the results of the performed evaluation on the bind join test cases with the 5 data stores and the 5 different selectivity factors. For each test case, the summary includes the execution times (in milliseconds) for performing traditional join (TJ) compared to bind join (BJ) and the performance improvement (PI = TJ/BJ). The test cases, for which there is no performance improvement, are omitted. Bind join produces an overhead, which is why for high values of the selectivity factor its usage is not advisable. First, when using bind join, the query engine cannot retrieve both tables in parallel, as it must wait for the outer table to be fully retrieved in order to get the set of distinct values of the outer keys. Second, it requires that this set of values be pushed into the subquery, which is handled differently by each data store, which explains the variability of performance improvements. The MonetDB wrapper makes a modification of the schema to create and populate a temporary table. MongoDB is doing faster, as it simply creates a filter out of the outer keys. With HBase, the process is much slower, as a native query is initiated by the query engine for each row of the outer table. The case of Sparksee is similar, as the Python code requests that a graph query be evaluated for each of the outer keys. And for ActivePivot, query rewriting takes place, however it brings a significant textual overhead, as the values of the outer keys must be specified with full dimension path.

D. Evaluation of Native Query Support

The support of native queries in CloudMdsQL allows the user to write powerful subqueries to the data stores that efficiently perform specific for the data store operations with native statements, rather than to express them in SQL and handle at the common query engine. This requires deeper expertise and knowledge about the specifics of the queried data stores; however, the added value is that the programmer can still request efficient operations at the data store, thanks to the support of native queries, and combine the results with other data store queries.

Let us consider the following queries, which request different operations to a data store, expressed in SQL. The queries focus on selection, aggregation, grouping, sorting and join operations. Let us assume that the programmer needs to request the same operations to non-SQL data stores (Sparksee and ActivePivot) that can natively support them. The naïve way to achieve this is for the wrapper developer or DBA to provide a static mapping of data store structures to relational tables and let the programmer express her query in SQL at the common query engine level. But now the programmer can take advantage of the native queries support and express the operations in native for the data store statements, which will result in an optimal execution. Note that all the generated TPC-H data for ActivePivot are stored in one single cube; so, although ActivePivot does not natively support joins between cubes, we categorize it as a join-capable data store only for this test case, because in fact a conceptual join between LINEITEM and ORDERS is handled by ActivePivot by dealing with dimensions of the same cube. With this group of test cases, we compare the execution times of the queries in two variants: when the requested operations are expressed in native queries and done by the data store (NQ), and when the operations are expressed in SQL queries and executed by the query engine after the wrapper simply delivers the entire table (SQ). Table 3 shows the comparison between NQ and SQ execution times, as well as the performance improvement (PI = SQ/NQ).

IV. CONCLUSION

This paper introduced the benchmarking of the CloudMdsQL polystore and evaluated the performance benefits in two specific categories of queries. First, the evaluation of bind join shows the ability of the query engine to run efficient queries across heterogeneous databases using both SQL and native approaches. Moreover, the measurements of the current or similar benchmark can be used as provision for a cost model based on inferring the performance improvement of bind joins for each particular data store. Second, the evaluation of the native query support stresses on the performance benefits provided by the high expressivity of the query language. With this benchmark we focused mostly on relational operators and their equivalents in non-SQL data stores. As a further work we will concentrate on evaluating the performance benefits in the context of data store specific operators, such as graph traversals or queries on nested documents.

  T2( ... WITHPARAMS orderkey long)@hbase = {* get 'orders', orderkey *} SELECT T1.*, T2(L_ORDERKEY).O_CUSTKEY, T2(L_ORDERKEY).O_ORDERSTATUS, ........

  CloudMdsQL allows named table expressions to be defined as Python functions, which is useful for querying data stores that have only API-based query interface. A Python expression yields tuples to its result set much like a user-defined table function. It can also use as input the result of other subqueries. Furthermore, named table expressions can be parameterized by declaring parameters in the expression's signature. For example, the following Python expression uses the intermediate data retrieved by T2 to return another table containing the number of occurrences of the parameter v in the array T2.z.

	T3(x int, c int WITHPARAMS v string)@python = {*
	for (x, z) in CloudMdsQL.T2:
	yield( x, z.count(v) )
	*}
	A (parameterized) named table can then be instantiated by
	passing actual parameter values from another native/Python
	expression, as a table function in a FROM clause, or even as a
	scalar function (e.g. in the SELECT list). Calling a named table
	as a scalar function is useful e.g. to express direct lookups into
	a key-value data store. In fact, in the current benchmark, we
	evaluate the usage of scalar lookups as an efficient alternative
	to the usage of expensive joins.

TABLE I

 I

	.	BIND JOIN SELECTIVITIES
	Filter predicate		Rows (approx.)	SF Rows / 1.5M
	L_PARTKEY = 10		24	0.002%
	L_PARTKEY <= 10		300	0.02%
	L_PARTKEY <= 100		3000	0.2%
	L_PARTKEY <= 1000		30000	2%
	L_PARTKEY <= 10000	300000	20%

TABLE II

 II

		.	BIND JOIN PERFORMANCE BENEFITS
	Data Store		0.002%	Selectivity factor 0.02% 0.2%	2%	20%
		TJ	56582	58113	59294	62027	67509
	MonetDB	BJ	953	1120	1982		8938	68757
		PI	59.4	51.9	29.9		6.9	1.0
		TJ	63783	65654	67807	75218	81146
	MongoDB	BJ	233	363	1124		6126	40908
		PI	273.7	180.9	60.3		12.3	2.0
		TJ	100576 102218 103340 105123
	HBase	BJ	1520	2711	10582	75886
		PI	66.2	37.7	9.8		1.4
		TJ	149694 151239 152116
	Sparksee	BJ	1480	7409	57739
		PI	101.1	20.4	2.6
		TJ	205071 206803 207944
	ActivePivot	BJ	1739	14409 145788
		PI	117.9	14.4	1.4

TABLE III .

 III NATIVE QUERY SUPPORT PERFORMANCE BENEFITS

	Data Store	Q_SEL	Q_AGG	Q_GRP	Q_ORD	Q_JOIN
		SQ	19068	26893	36870	51075	69393
	Sparksee	NQ	76	326	2654	2640	785
		PI	250.9	82.5	13.9	19.3	88.4
	Active Pivot	SQ NQ PI	20009 261 76.7	12592 112 112.4	10067 16 629.2	10907 77 141.6	193963 189 1026.3

This research has been partially funded by the European Commission under projects CoherentPaaS, LeanBigData, and CloudDBAppliance (grants FP7-611068, FP7-619606, and H2020-732051), the Madrid Regional Council, FSE and FEDER, project Cloud4BigData (grant S2013TIC-2894), the Spanish Research Agency MICIN project BigDataPaaS (grant TIN2013-46883), and the Spanish CDTI/MEC NEOTEC program (grant SNEO-20151285).