[
DG method for acoustoelastic waves ]
A Discontinuous Galerkin method with a modified
penalty flux for the propagation and
scattering of acoustoelastic waves
R. Ye et al.]
Ruichao Ye, Maarten V. de Hoop, Christopher L. Petrovitch,
Laura J. PyrakNolte, Lucas C. Wilcox
Department of Earth Science, Rice University, Houston TX, USA
Simons Chair in Computational and Applied Mathematics and Earth Science, Rice University, Houston TX, USA
Applied Research Associates, Inc., RaleighDurham NC, USA
Department of Physics, Purdue University, West Lafayette IN, USA
Department of Applied Mathematics,
Naval Postgraduate School, Monterey CA, USA
We develop an approach for simulating acoustoelastic wave phenomena, including scattering from fluidsolid boundaries, where the solid is allowed to be anisotropic, with the Discontinuous Galerkin method. We use a coupled firstorder elastic strainvelocity, acoustic velocitypressure formulation, and append penalty terms based on interior boundary continuity conditions to the numerical (central) flux so that the consistency condition holds for the discretized Discontinuous Galerkin weak formulation. We incorporate the fluidsolid boundaries through these penalty terms and obtain a stable algorithm. Our approach avoids the diagonalization into polarized wave constituents such as in the approach based on solving elementwise Riemann problems.
iscontinuous Galerkin method – penalty flux – fluidsolid boundaries – anisotropy
1 Introduction
The accurate computation of waves in realistic threedimensional Earth models represents an ongoing challenge in local, regional, and global seismology. Here, we focus on simulating coupled acoustoelastic wave phenomena including scattering from fluidsolid boundaries, where the solid is allowed to be anisotropic, with the Discontinuous Galerkin method. Of particular interest are applications in geophysics, namely, marine seismic exploration and global Earth inverse problems using earthquakegenerated seismic waves as the probing field. In the first application, we are concerned with the presence of the ocean bottom and in the second one with the coremantleboundary (CMB) and innercoreboundary (ICB). Our formulation closely follows the analysis of existence of (weak) solutions of hyperbolic firstorder systems of equations by [Blazek et al.(2013)Blazek, Stolk, & Symes]. We use an unstructured tetrahedral mesh with local refinement to accommodate highly heterogeneous media and complex geometries, which is also an underlying motivation for employing the Discontinuous Galerkin method from a computational point of view.
In the past three decades, a wide variety of numerical techniques has been employed in the development of computational methods for simulating seismic waves. The most widely used one is based on the finite difference method [e.g., [Madariaga(1976)] and [Virieux(1986)]]. This method has been applied to computing the wavefield in threedimensional local and regional models [e.g., [Graves(1996)] and [Ohminato & Chouet(1997)]]. The use of optimal or compact finitedifference operators has provided a certain improvement [e.g., [Zingg et al.(1996)Zingg, Lomax, & Jurgens] and [Zingg(2000)]]. Methods that resort to spectral and pseudospectral techniques based on global gridding of the model have also been used both in regional [e.g., [Carcione(1994)]] and global [e.g., [Tessmer et al.(1992)Tessmer, Kessler, Kosloff, & Behle] and [Igel(1999)]] seismic wave propagation and scattering problems. However, because of the use of global basis functions (polynomial: Chebyshev or Legendre, or harmonic: Fourier), these techniques are limited to coefficients which are (piecewise) sufficiently smooth. The finite difference method suffers from a limited accuracy in the presence of a free surface or surface discontinuities with topography within the model [e.g., [Robertsson(1996)] and [Symes & Vdovina(2009)]]. A procedure for the stable imposition of freesurface boundary conditions for a secondorder formulation can be found in [Appelö & Petersson(2009)]. Another approach, belonging to a broader family of interface methods, handles both free surfaces [e.g., [Lombard et al.(2008)Lombard, Piraux, Gélis, & Virieux]] and fluidsolid interfaces [e.g., [Lombard & Piraux(2004)]] in such a way, conjectured by the authors, that enables higherorder accuracy to be obtained. [Kozdon et al.(2013)Kozdon, Dunham, & Nordstræm] use summationbyparts finite difference operators along with a weak enforcement of boundary conditions to develop a multiblock finite difference scheme which achieves higherorder accuracy for complex geometries.
A key development in the computation of seismic waves has been based on the spectral element method (SEM) [[Komatitsch & Tromp(2002)]]. In its original formulation, in terms of displacement [[Komatitsch & Vilotte(1998)]], continuity of displacement and velocity is enforced everywhere within the model. In the case of a boundary between an inviscid fluid and a solid, however, the kinematic boundary condition is perfect slip; therefore, only the normal component of velocity is continuous across such a boundary, and thus this formulation is not applicable. Some classical finiteelement methods (FEMs) alternatively introduce coupling conditions on fluidsolid interfaces between displacement in the solid and pressure in the fluid [e.g. [Zienkiewicz & Bettess(1978), Bermúdez et al.(1999)Bermúdez, HervellaNieto, & Rodrıguez]].
The FEM and SEM are commonly (but not exclusively) based on the secondorder form of the system of equations describing acoustoelastic waves. In this case, the acoustoelastic interaction is effected by coupling the respective wave equations through appropriate interface conditions. To resolve the coupling, a predictormulticorrector iteration at each time step has been used [[Komatitsch et al.(2000)Komatitsch, Barnes, & Tromp], [Chaljub et al.(2003)Chaljub, Capdeville, & Vilotte]]. A computationally more efficient time stepping method for global seismic wave propagation accommodating the effects of fluidsolid boundaries, as well as transverse isotropy with a radial symmetry axis and radial models of attenuation, was proposed in [Komatitsch et al.(2005)Komatitsch, Tsuboi, & Tromp]. It uses a velocity potential formulation and a secondorder accurate Newmark time integration, in which a time step is first performed in the acoustic fluid and then in the elastic solid using interface values based on the fluid solution. Currently the SEM is used in a variety of implementations in global and regional seismic simulation, with the effects of variations in elastic parameters, density, ellipticity, topography and bathymetry, fluidsolid interfaces, anisotropy, and selfgravitation included [e.g. [Carrington et al.(2008)Carrington, Komatitsch, Laurenzano, Tikir, Michéa, Le Goff, Snavely, & Tromp]].
In contrast to classical finite element discretizations, the Discontinous Galerkin (DG) method imposes continuity of approximate solutions between elements only weakly through a numerical flux. The Discontinuous Galerkin method has been employed for solving secondorder wave equations in both the acoustic and elastodynamic settings [e.g. [Rivière & Wheeler(2000)], [Grote et al.(2006)Grote, Schneebeli, & Schötzau], [Chung & Engquist(2006)] and [De Basabe et al.(2008)De Basabe, Sen, & Wheeler]]. [Etienne et al.(2010)Etienne, Chaljub, Virieux, & Glinsky] employ a central numerical flux in a DG scheme combined with a leapfrog time integration for the velocitystress elasticwave formulation. [Dumbser & Käser(2006), Käser & Dumbser(2008)] developed a nonconservative formulation with an upwind numerical flux using only the material properties from the side of the interface that is opposite to the outer normal direction. [Wilcox et al.(2010)Wilcox, Stadler, Burstedde, & Ghattas] derived an upwind numerical flux by solving the exact Riemann problem on interior boundaries of each element with material discontinuities based on a velocitystrain formulation of the coupled acoustoelastic equations.
In this work, we essentially extend the upwind flux, given by [Warburton(2013)] for hyperbolic systems, to a penalty flux based on the boundary continuity condition for general fluidsolid interfaces. The novelties of our approach are the following: we

use a coupled firstorder elastic strainvelocity, acoustic velocitypressure formulation,

obtain a selfconsistent Discontinuous Galerkin weak formulation without diagonalization into polarized wave constituents,

append penalty terms, derived from interior boundary continuity conditions, with an appropriate weight to the numerical (central) flux so that the consistency condition holds for the discretized Discontinuous Galerkin weak formulation,

incorporate fluidsolid boundaries through the mentioned penalty terms.
We note that the DG method is naturally adapted to wellposedness, in the sense that it makes use of coercivity of the operator defining the part of the system containing the spatial derivatives separately in the solid and fluid regions.
2 The system of equations describing acoustoelastic waves
We consider a bounded domain which is divided into solid and fluid regions, and , respectively. The interior boundaries include solidsolid interface , fluidfluid interface , and fluidsolid interface , (where we distinguish whether the fluid or solid is on a particular side). We present the weak form of the coupled acoustoelastic system of equations.
Hooke’s law in an elastodynamical system is expressed by relating stress, , and strain, . Assuming small deformations gives a linear relationship, that is, , where is the stiffness tensor. Through the relevant symmetries, this tensor only contains 21 independent components. We use the Voigt notation which simplifies the writing of tensors while introducing and . In this notation the stiffness tensor takes the form of a 6 by 6 matrix, C, defined by,
(1) 
The isotropic case is obtained by setting all of the components to zero except for , , , , and ; are the Lamé parameters. Furthermore, denotes the density. The anisotropic elastodynamical equations are written in terms of the strain, E, and the particle velocity, v,
(2) 
in . In fluid regions, , we use the pressurevelocity formulation,
(3) 
Here, is the pressure, while we use to distinguish acoustic field quantities and material parameters from the elastic ones. In the above, denotes a volume source density of injection and f denotes a volume source density of force.
The solidsolid, fluidsolid and fluidfluid boundary conditions are given by
(4a)  
(4b)  
(4c) 
The convention is determined by the direction of the interface normal, n. The outer normal vector points in the direction of the “” side of the interface.
We introduce test functions (tensors) in the solid regions and in the fluid regions, which are assumed to be contained in the same spaces and satisfy the same boundary conditions as and . Using (2) and (3), we find that
(5a)  
(5b)  
(5c)  
(5d) 
Assuming an outer tractionfree boundary condition in (5b) and an outer pressurefree boundary condition in (5c), and applying an integration by parts, we obtain
(6a)  
(6b) 
We use the fluidsolid boundary conditions (4b), replacing the fluidsolid surface integrals in (6a) and (6b) by taking the average of both sides consistent with a central flux scheme, and obtain
(7a)  
(7b) 
This form of the equations is analogous to the one used in the spectral element method, see [Chaljub & Valette(2004)]. Applying an integration by parts, again, in (7), we recover the coupled strong formulation,
(8a)  
(8b) 
We use this system of equations together with (5a) and (5d) to develop our Discontinuous Galerkin method based approach.
3 Discontinuous Galerkin method with fluidsolid boundaries
The domain is partitioned into elements, . We distinguish elements, , in the solid regions from elements, , in the fluid regions. Correspondingly, we distinguish fluidfluid (), solidsolid () and fluidsolid () faces for each element; thus the interior boundaries are decomposed as
and so are the elements’ boundaries: and . The mesh size, , is defined as the maximum radius of each tetrahedral’s inscribed sphere.
We introduce the broken polynomial space where the local space is defined elementwise as , with a set of polynomial basis further discussed in Section 3.2. The subscript “” indicates the refinement of with decrease in mesh size. The semidiscrete timedomain, discontinuous Galerkin formulation using a central flux yields: Find , with each component for each one of them in such that
(9) 
and
(10) 
hold for each element or , for all test functions . The notations and indicate polynomial approximation of f and . Here,
(11a)  
(11b) 
in the solid regions, while
(12a)  
(12b) 
in the fluid regions, using interior boundary continuity conditions. A similar formulation for Maxwell’s equations, using the central flux, can be found in [Hesthaven & Warburton(2007), Chapter 10, Page 434].
3.1 Energy function of central flux
We consider a timedependent energy function comprising both the solid and fluid regions, , with
(13) 
The functions in (13) define a norm both in the solid and in the fluid regions. Taking the time derivative and noting that C is symmetric, we have
(14)  
(15) 
Starting from (9) and (10) and carrying out the summation over all the elements yields
(16) 
This property is obtained as follows:
In (9) and (10) we let , and obtain elementwise
(17) 
and similarily
(18) 
From (9), (11), (14) and (17),
()  
() 
In the above,
(19) 
The surface integration terms cancel out when summed from both sides of the solidsolid interfaces because of the continuity condition (4a) and the opposite outer normal directions. We are left with the contributions from solidfluid inner faces, ,
(20) 
A similar result in the fluid region obtained from (10), (12), (15) and (18) yields
(21) 
and the surface integration terms on the solidfluid and fluidsolid interfaces in (20) and (21) cancel out due to (4b). Therefore (16) is obtained. We note that the surface integration along solidfluid interfaces and are essential to guarantee energy conservation.
3.2 Nodal basis functions
The discretized solution follows an expansion, componentwise, into nodal trial basis functions of order , as is in [Hesthaven & Warburton(2007)],
(22) 
and similarly for the other fields, . The superscript, , indicates a local expansion within element . In the above, is a set of threedimensional Lagrange polynomials associated with the nodal points, (see Figure 1), with each polynomial defined as
We use the warp & blend method [[Warburton(2006)]] to determine the coordinates of nodal points in the tetrahedron by numerically minimizing the Lebesgue constant of interpolation. For an order interpolation there are nodal points.
= 1  = 3  = 8 
The medium coefficients are expanded in a likewise manner
(23) 
and similarly for . When refining a mesh, we expect an increase in number of elements with decreased size.
3.3 The system of equations in matrix form
To simplify the notation in the further development of a numerical scheme, we introduce a joint matrix form of the system of equations. We map the components of and to and matrices, respectively,
(24) 
and, correspondingly, the components of body forces f and to the matrix
Equations (2) and (3) attain the form
(25) 
where
and
that is,
with