Explicit Polynomial Trefftz-DG Method for Space-Time Elasto-Acoustics

The popularity of Discontinuous Galerkin (DG) methods has increased during the last years in particular with easy access to high performance computers. DG approximations gather different interesting features such as high-order accuracy, hp-adaptivity, easy parallelization but have the disadvantage of requiring a high number of degrees of freedom. DG methods are thus known to use higher computational burden and herein, we investigate the possibility of reducing it by employing an appropriate set of basis functions. We then solve the first-order elasto-acoustic wave problem in a space-time Trefftz-DG space. We use a space-time mesh with space-time polynomial basis functions. The space-time variational formulation reduces then to surface integrals set only on the skeleton of the mesh. However, the solution of the problem requires inverting a sparse matrix whose size could be very large when addressing realistic geophysical problems. A direct implementation of the Trefftz-DG formulation creates thus extra cost as compared to that of standard DG formulation where the space and time integrations are performed one by one after another. We then propose an approximate formulation which amounts to inverting the block-diagonal component of the Trefftz-DG matrix. We then end up with an explicit representation of the wavefield which is stable provided a CFL condition is satisfied.

In this section, following [10] and the framework therein, we propose a formulation of the elasto-acoustic coupling reading as a first-order system. Here and further the sub-scripts F and S corresponds to the acoustic (fluid) and elastodynamic (solid) domains.

Elasto-Acoustic Equations
We introduce a space-time domain Q ≡ ( F ∪ S ) × I , where F ⊂ R d is a bounded Lipschitz domain of dimension d filled with fluid, S ⊂ R d is a bounded Lipschitz elastodynamic domain of dimension d filled with solid, and I ≡ [0, T ] is the time interval. All medium parameters c F ≡ c F (x) and ρ F ≡ ρ F (x), standing for the acoustic wave propagation velocity and fluid density respectively, as well as the inverted stiffness tensor C −1 (x) ≡ A(x) and the solid density ρ S ≡ ρ S (x), are assumed to be piecewise constant and positive. We denote by F S = F ∩ S the fluid-solid interface. The elasto-acoustic system of equations is based on the coupling of the first-order acoustic equation, written in terms of velocity v F ≡ v F (x, t) and pressure p ≡ p(x, t) fields: where n F is the normal vector to ∂ F , the source term f ≡ f (x, t), the boundary condition g F , the velocity v F 0 and the pressure p 0 are the initial data, with the firstorder elastodynamic system, written in terms of velocity v S ≡ v S (x, t) and stress tensor (symmetrical and positive) σ ≡ σ (x, t) fields: where n S is the normal vector to ∂ S , the boundary condition g F , the velocity v S 0 and the stress tensor σ 0 are the initial data. The transmission conditions between the two systems (2) and (1) represent the continuity of velocity and stress normal components F S : The velocities aligned with the interface and the tangential stress remain unconstrained.

Space-Time Trefftz-DG Formulation
We introduce a non-overlapping space-time mesh T h on Q composed of space-time Lipschitz elements K F ⊂ F × I and K S ⊂ S × I . We denote by T F h (resp. T Sh ) the restriction of T h to the fluid (resp. solid) domain. Let n K F ≡ (n x K F , n t K F ) be the outward-pointing unit normal vector on ∂K F , and n K S ≡ (n x K S , n t K S ) be the outward-pointing unit normal vector on ∂K S . We assume that all medium parameters are constant in K F and K S respectively. The mesh skeleton F h ≡ We consider the test functions ω F , q, ω S , ξ in for v F , p, v S and σ respectively, where the Trefftz space T(T h ) is defined on the mesh T h as follows: This space is of Trefftz type since it is a subspace of the regular space V h (T h ) composed of local solutions of the volumic governing equations (1) and (2) set in each element K F and K S respectively. As in the standard DG methods, the next step in order to obtain the variational formulation consists in multiplying the equations of (1) by the test functions q and ω F in T(T h ), and the equations of (2) by the test functions ξ and ω S in T(T h ) respectively, and, as is standard in space-time DG methods, we integrate by parts the obtained equations not only in space but also in time: Thanks to the choice of test functions the left hand side of the above space-time formulation contains only surface integrals. The numerical fluxes in timev F h ,p h , v Sh ,σ h and in spacev F h ,p h ,v Sh ,σ h are defined in the standard DG notations [2,3,7] as follows: ⎛ Here, α 1 , α 2 , β 1 , β 2 , δ 1 , δ 2 , γ 1 , and γ 2 are positive penalty parameters. As in standard DG methods, a suitable choice of these penalty parameters allows one to prove stability of the overall method. It is shown in [2,3] that they contribute to the accuracy and convergence of the numerical method. We refer to [2,3] for more details on the definition of the numerical fluxes.
Summing the contribution (4) of all elements K F , K S ∈ T h , and introducing the bilinear A T DG (· ; ·) and the linear T DG ·) forms for the left-hand side and the right-hand side expressions respectively, we obtain the Trefftz-DG formulation for the elasto-acoustic problem: The analysis of well-posedness of (5) is based on the coercivity and continuity estimates of the bilinear and linear forms in mesh-dependent norms [2,3]. The proof is similar to the one given in [10] where the acoustic wave equation is addressed. In Sect. 2 we provide the algorithm of the Trefftz-DG formulation (5), and we discuss different analytical and numerical approaches for its optimization.

Implementation of the Algorithm
The numerical implementation of the Trefftz-DG formulation is different from the standard DG ones which address the space and time integration separately. Standard DG space integrations have the interesting feature of leading to a block-diagonal mass matrix and allow then the use of explicit time integration. The computational costs thus depend on a CFL condition which sets the value of the time step as a function of the space step. On the other hand, a naive implementation of Trefftz-DG methods require performing a space-time integration which leads to invert a sparse matrix whose size tends to be huge. It is thus not obvious that a crude implementation of the Trefftz-DG algorithm does not generate additional cost as compared to standard DG ones.
In this section we provide some important steps of implementation of Trefftz-DG formulation (5) and discuss optimization techniques. The complete algorithm with more numerical details can be found in [2,3].

Change-Over Between the Time Slabs
To simplify the presentation, we assume here that we use the same order of approximation on each cells, so that we have N When compared to the computational cost of standard DG implementation, the corresponding Trefftz-DG cost is thus increased and it is mainly due to the need of inverting the large-sized matrix. The most obvious way to reduce the size of the matrix, which is classically used in most work on space-time Trefftz method, is to consider time slabs. We restrict ourselves to the case of cartesian meshes, but this methodology can also be applied to unstructured meshed. An alternative is to use tent-pitched meshes that respect the causality, this will be the topic of a future work. In order to optimize the execution of the algorithm, we propose to divide the space-time domain Q into N t elementary time slabs Q 1 , Q 2 , . . . , Q N t and to solve the problem slab by slab, considering the final results, computed in the current time slab at time t, as initial values for the next slab at time t + t (see Fig. 2). Thus the size of matrix inside each time slab is N t times smaller, compared to the initial one. Moreover, if the medium parameters are fixed in time, and the space discretization is preserved from slab to slab, the matrix can be computed and inverted once, and then re-projected onto the next time slabs, reducing thus the global numerical cost.

Polynomial Basis
One of the important advantages of Trefftz type methods is the flexibility in the choice of basis functions provided they satisfy the Trefftz property locally in each element. To perform the numerical simulations, we have extended the algorithm proposed by Maciag in [9] for computing wave polynomials, solutions of the second order transient wave equation, to the first order acoustic and elastodynamic systems of dimension one and higher. It consists in computing a polynomial basis, defined in the reference element, using Taylor expansions of generating exponential functions which are local solutions of the initial system of equations. An example of spacetime wave polynomial basis for the first-order acoustic wave equation reads as follows (approximation degree p=3, dimension of the physical space d = 1): This basis contains the couples of polynomial functions (φ v · ,φ p · ), corresponding to the velocity and pressure respectively, which are locally defined and satisfy the Trefftz property inside each element of the mesh, and of degrees less or equal to p (p = 0, 1, 2, 3) to provide an approximation of order p. By their construction, the Trefftz basis functions are not attached to the coordinates of the degrees of freedom inside the element, contrary to the Lagrange polynomials. Even if we compute only surface integrals, we can evaluate the final approximation solution in any point of the element refinement. We refer to [2] for more numerical details as well as for the acoustic and elastodynamic basis examples of higher dimensions.

Inversion of the Matrix M Inside a Time Slab
The inversion of the matrix inside the time slab can be explicitly reduced to the inversion of its block-diagonal component, which corresponds to the integration at the bottom and top of the time slab (initial and final time faces F 0 h and F T h ), thanks to the Taylor expansion formulas. More precisely, let us recall the expression for the bilinear form A T DG (· ; ·) from Sect. 1.2: Table 1 Accuracy (L 2 -error in space and in time) of the solution when using the approximate inversion with n = 3, 4, 5 and the exact inversion The accelerating factor is the ratio of the computational costs of the two methods for reaching the same accuracy

Numerical Tests
For the numerical implementation of the Trefftz-DG method we have considered a 2D medium composed of two homogeneous rectangular layers: the acoustic one and the elastodynamic one. We have set a source term at the fluid-solid interface, and two receivers in the acoustic layer and in the elastodynamic one. The numerical signals at both receivers have been validated with the analytical solutions computed with Gar6more code [5]. In Fig. 3 we show the convergence of the numerical velocity as a function of cell size for different degrees of approximation (p = 0, 1, 2, 3) computed at receivers in (a) 2D acoustic layer and (b) 2D elastodynamic layer. In each case, the convergence rate is higher than the corresponding approximation degree. We refer to [2], where we provide more examples.

Conclusion
The Trefftz-DG methodology for solving the first order elasto-acoustic system has demonstrated the important advantages, such as the use of degrees of freedom evaluated at the element faces only, the flexibility in the choice of the basis functions and the unconditional stability. However, in its initial form, it still shows some limitations due to the space-time integration that leads to the representation of the discrete system by a huge sparse matrix whose straightforward inversion is very expensive, even when using time slabs. We find ourselves in a situation of using an implicit scheme for solving the forward problem that risks to overload the iterative process of the corresponding inverse problem in order to reconstruct very large propagation domains. Fortunately, thanks to the decomposition of the matrix by separating the time variables from the space ones, we could benefit from the block-diagonal structure of the standard DG formulation ending up with an explicit scheme, that is more convenient from the numerical point of view. The performed numerical tests clearly illustrate the interest of the split version of discrete problem. The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.