Simulations of wave propagation through elastic-acoustic coupling media are applicable to a wide range of scientific and engineering areas. For example, coupled elastic-acoustic media arises when simulating wave propagation through the human bone and tissue. While wave propagation in tissue is modeled by the acoustic wave equation, wave propagation in bone is more accurately modeled using the elastic wave equation, and when considering wave propagation through both bone and tissue, careful attention is required for treatment of the elastic-acoustic interface. Wave propagation through coupled elastic-acoustic media also arises in seismology, where oceans are modeled as acoustic materials and the earth is modeled as an elastic medium.
Several high order finite element methods have been developed for coupled elastic-acoustic wave propagation based on both first and second order formulations of the underlying equations. In [4], Komatitsch et al. use a spectral element method (SEM) for the second order form of the equations, and enforce the coupling between acoustic and elastic media using with a predictor-multicorrector iteration at each time step. A more efficient time-stepping approach based on explicit coupling conditions is proposed in [5,6]. Discontinuous Galerkin (DG) methods have also been developed for coupled elasticacoustic media, with elastic-acoustic interface conditions typically incorporated through modifications of the numerical flux. For second order equations, Antonietti et al. [7] analyze the stability and convergence of a symmetric interior penalty DG formulation on polygonal and polyhedral meshes. Appelo and Wang [8]i n t r o d u c e an "energy-based" second order DG method which can be made to either conserve or dissipate energy based on the choice of numerical flux.
DG methods, which were originally developed for first-order hyperbolic equations, are also widely used for first-order formulations. Wilcox et al. [9]d e r i v e an upwind numerical flux from the exact Riemann problem at acoustic-acoustic, elastic-elastic, and elastic-acoustic interfaces, and use this to construct a first-order velocity-strain DG-SEM scheme for coupled isotropic elastic-acoustic media on meshes of curved hexahedral elements. The authors show stability of the continuous DG formulation; however, semi-discrete stability in the presence of inexact quadrature, curved meshes, and sub-cell heterogeneities is not discussed in detail. In [10], Zhan et al. extend this DG-SEM method to anisotropic elastic-acoustic media by solving a simplified Riemann problem on inter-element interfaces, though high order accuracy and energy stability are not addressed theoretically. Ye et al. [11]c i r c u m v e n t the Riemann problem altogether by using a DG formulation with a dissipative upwind-like "penalty" flux. The resulting DG method is high order accurate and provably energy stable for anisotropic elastic-acoustic media with piecewise constant heterogeneities.
In this paper, we develop a high order DG method for coupled elastic-acoustic media based on the first-order stressvelocity form of the equations. The proposed method utilizes a simple dissipative upwind-like penalty flux and weightadjusted mass matrices (a generalization of mass lumping) [1][2][3]. The method applicable to unstructured and curved tetrahedral meshes, and is high order accurate and energy stable in the presence of arbitrary heterogeneous media including anisotropy and micro (sub-cell) heterogeneities. Instead of an exact upwind flux, we add upwind-like dissipation through a penalty flux based on natural continuity conditions between acoustic-acoustic, elastic-elastic, and coupled elasticacoustic interfaces. Like the upwind flux, the penalty flux adds dissipation and achieves theoretically optimal high order convergence rates for all numerical experiments without impacting the maximum time-step size. However, expressions for the penalty flux are significantly simpler than the fluxes developed by Wilcox et al. and Zhan et al. [9,10]. Additionally, we prove that the penalty flux is consistent and that the semi-discrete DG formulation is energy stable for general "modal" DG formulations in the presence of both sub-cell heterogeneities and curved elements. Experiments with high order DG discretizations on curvilinear simplicial meshes verify these theoretical properties.
The outline of the paper is as follows: In Section 2, we review DG formulations for the acoustic and elastic wave equations. In Section 3, we introduce the numerical flux for elastic-acoustic interfaces and prove that the resulting DG formulation is energy stable and consistent. In Section 4, we verify the stability and accuracy of the proposed DG method, and conclude in Section 5 with an application in photoacoustic tomography (PAT).
In this section, we briefly review high order DG discretizations for the acoustic and elastic wave equations. In the presence of micro (sub-cell) heterogeneities, inverse weighted mass matrices appear in the matrix forms of these discretizations. These inverses are approximated using easily invertible weight-adjusted mass matrices, resulting in a weight-adjusted DG method. The weight-adjusted approach will be extended to the DG formulation for coupled elastic-acoustic wave propagation and curvilinear meshes in Sections 3.1 and 3.3.
We assume a physical domain , which is exactly represented by a triangulation h consisting of K non-overlapping elements D k . We assume that each element D k is the image of the reference element D under a mapping k
where x = (x, y, z) are physical coordinates on the kth element and x = ( x, y, z) are coordinates on the reference element.
Over each element D k , we define the polynomial approximation space V h D k as
where V h D is a polynomial approximation space of degree N on the reference element. In this work, 1 the reference element is taken to be bi-unit right triangle,
and the reference approximation space V h D is taken to be total degree N polynomials,
1 In three dimensions, the reference element and approximation space are the bi-unit right tetrahedron and total degree N polynomials
On an element D k , we define the jump of scalar and vector valued functions across element interfaces as
where p + , u + and p, u are the neighboring and local traces of the solution over the interface, respectively. Note that, for a shared interface between two elements D k and D k,+ , the sign of the jump is different depending on whether the jump is defined with respect to D k or D k,+ . The average across an interface is defined as
In this work, we use a first-order pressure-velocity formulation for the acoustic wave equation (e.g. in fluid media)
where p is the acoustic pressure, u ∈ R d is the vector of velocities in each coordinate direction, and ρ and c are density and wavespeed, respectively. For simplicity, we assume unit density ρ = 1. We also assume that (1)i s posed over time t ∈[0, T )
on the physical domain with boundary ∂ , with the wavespeed bounded from above and below by
We adopt a DG variational formulation from [12], which is given over element D k by 1
where n is the outward normal vector, and τ p , τ u are penalty parameters. Here, (u, v) L 2 D k and u, v L 2 ( f ) denote the L 2 inner products over D k and a face f of the surface ∂ D k , respectively.
For the elastic wave equation, we use a symmetrized first-order stress-velocity formulation from [3]. Let ρ be the density and C be the symmetric matrix form of constitutive tensor relating stress and strain. The first-order system in d dimensions is given by
where v is the vector of velocity and σ is a vector consisting of unique entries of the symmetric stress tensor. In two dimensions, the matrices A i are defined as
and the expression of A i in three dimensions can be found in [3]. Note that, by factoring out C , the resulting matrices A i do not involve any material coefficients and all entries are either 0o r 1. The elastic wave equation is discretized using the following DG formulation:
where A n is normal matrix defined as A n = d i=1 n i A i , and terms τ v , τ σ are penalty parameters introduced on element interfaces. The DG formulations (2) and (4)a r e provably consistent and energy stable for non-negative penalty parameters τ p , τ u , τ v , τ σ ≥ 0 [ 1,3].
The matrix form of the DG formulations in the previous section involve mass and differentiation matrices. We assume the reference and physical approximation spaces V h D and V h D k are spanned by bases {φ i } N p i=1 and {φ k i } N p i=1 , respectively. The mass matrix M k , weighted mass matrix M k w and face mass matrix M k f for D k are defined as
where J k and J k f are the volume and face Jacobian of the affine mapping k , and w(x) is a spatially varying positive and bounded weight. We also define weak differentiation matrices S i with entries
Using the above notation, the DG formulation (2)c a n be written in matrix form as
where U i and p are degrees of freedom for u i and p. The flux terms F p , F u are defined such that
The DG scheme (4)f o r the elastic wave equations can similarly be written as
where F v , F σ denote the elastic flux terms, ⊗ denotes the Kronecker product, and the matrix-valued weight mass matrix
denotes the scalar weighted mass matrix with weight C -1 ij .
In this work, we pair high order DG methods with explicit time-stepping schemes, which require the inversion of DG mass matrices at each time-step. Let U denote the vector of all DG degrees of freedom, and let M k 1/c 2 , A k denote the local matrices representing the local DG mass matrix and spatial DG formulation, such that the semi-discrete DG scheme can be written over D k as follows:
When the wavespeed c 2 is approximated by a constant over each element, it is possible to apply M k 1/c 2 -1 using only the constant values of J k , c 2 over each element and a single reference mass matrix inverse M -1 over the entire mesh. However, inverses of weighted mass matrices are distinct from element to element when c 2 possesses sub-element variations. Typical implementations precompute and store these weighted mass matrix inverses [13,14], which significantly increases the storage cost of high order DG schemes.
To address this issue, we use a weight-adjusted discontinuous Galerkin (WADG) proposed in [1,3], which is energy stable and high order accurate for sufficiently regular weighting functions. WADG approximates each weighted mass matrix by a weight-adjusted approximation
Since the weight only appears in
can be applied using reference inverse mass matrices and a matrix-free quadrature-based evaluation of M k 1/w . Analogously, the inverse of M C -1 can be approximated by the inverse of a matrixweighted weight-adjusted mass matrix
In practice, weight-adjusted mass matrix inverses are applied in a matrix-free fashion using sufficiently accurate quadrature rules. We follow [1] and use simplicial quadratures which are exact for polynomials of degree 2N + 1 [ 15]. Let x i , w i denote the quadrature points and weights on the reference element D . We define the interpolation matrix V q as
whose columns consist of values of basis functions at quadrature points. On each element D k , we have
where k x i are quadrature points on D k and c 2 k x denote the values of the wavespeed at quadrature points. Plugging the approximation (6)i n t o the local DG formulation (5), we obtain
Evaluating M k -1 A k U is equivalent to the evaluation of the DG right hand side for a unit weight 1/c 2 = 1. Evaluating the remainder of the right hand side of (7)r e q u i r e s applying the product of an unweighted mass matrix and weighted mass matrix. This can be done using quadrature-based matrices as follows:
where P q = M -1 V T q diag ( w) is a quadrature discretization of the polynomial L 2 projection operator on the reference el- ement. Moreover, since P q , V q are reference operators, the implementation of (8)r e q u i r e s only O N d storage for values of the wavespeed c 2 k x at quadrature points for each element. In contrast, storing full weighted mass matrix inverses or factorizations requires O N 2d storage on each element. For example, in three dimensions, the number of quadrature points on one element scales with O (N p ) = O (N 3 ), while number of entries in each weighted mass matrix inverse is
For the first-order acoustic and elastic wave equations, the discontinuous Galerkin schemes ( 2) and ( 4)a r e consistent and discretely energy stable for a large class of quadrature rules. The goal of this work is to extend these existing schemes to solve wave problems in coupled elastic-acoustic media. The challenge is to derive an appropriate numerical flux for the interface between acoustic and elastic domains. In this section, we propose a new numerical flux across elastic-acoustic interfaces, and prove the consistency and discrete energy stability of the elastic-acoustic DG formulation under this new flux.
We begin with the continuity conditions on the interface between different media. For an acoustic-acoustic interface, the normal velocity and pressure are continuous, i.e.,
For an elastic-elastic interface, the velocity and the traction are continuous, i.e.,
For an interface between elastic and acoustic media, the normal component of the velocity and the traction are continuous, i.e.,
where u and v denote velocity in acoustic and elastic media, respectively. Based on these continuity conditions, we derive an upwind-like numerical flux for the elastic-acoustic interface.
For clarity, we will distinguish between acoustic and elastic fluxes at a coupled elastic-acoustic interface. Let e h , a h denote the elastic and acoustic computational domains, respectively. Let Ŵ ea and Ŵ ae denote the respective boundaries of e h and e h which correspond to the elastic-acoustic interface. Let Ŵ aa , Ŵ ee be the complement of Ŵ ae and Ŵ ea in collections of all acoustic element faces and elastic element faces, respectively. On Ŵ ae , the numerical fluxes are taken to be 1 2
while the numerical fluxes on Ŵ ea are given by 1 2
We now formulate a DG scheme for the first-order elastic-acoustic coupled wave equations. In the acoustic domain a h , the DG formulation is given by 1
In the elastic domain e h , the DG formulation is given by
We note that media heterogeneities are incorporated into the left hand side of the DG formulations ( 10) and ( 11), and that the numerical fluxes are independent of any variations in 1/cfoot_0 , C -1 . In our numerical experiments, we approximate the weighted mass matrices induced by micro (sub-cell) heterogeneities in 1/c 2 , C -1 by easily invertible weight-adjusted mass matrices as described in Section 2.4.
In this section, we prove that the DG formulations ( 10) and ( 11)a r e consistent and energy stable in arbitrary heterogeneous media.
Theorem 3.1. The coupled discontinuous Galerkin scheme is consistent.
Proof. Assume that u, p, v, σ are exact solutions of coupled elastic-acoustic wave equations, and that boundary conditions are imposed through consistent modifications of the numerical flux. 2 Then, plugging them into (10) and ( 11)c a u s e s the volume terms to vanish. Consistency follows if the numerical flux terms also vanish.
At acoustic-acoustic interfaces, the pressure and normal velocity are continuous. Thus, the numerical flux reduces to 1 2
At elastic-elastic interfaces, the traction A T n σ and the velocity are continuous, and the numerical flux reduces to
For an elastic-acoustic interface Ŵ ae , we have
Similarly, on Ŵ ea , we have
Thus, consistency holds for acoustic-acoustic, elastic-elastic and elastic-acoustic interfaces, which implies the coupled DG scheme is consistent.
The formulations (10) and ( 11)c a n also be shown to be energy stable for any choice of τ u = τ v ≥ 0, τ p = τ σ ≥ 0. For simplicity, we assume zero homogeneous Dirichlet boundary conditions on ∂ in the proof of energy stability.
Theorem 3.2. The coupled discontinuous Galerkin scheme is energy stable for τ u = τ v ≥ 0, τ p = τ σ ≥ 0, in the sense that
where a h and e h denote the acoustic and elastic computational domain, respectively.
Proof. For the acoustic part, taking q = p, w = u and integrating the divergence term of the pressure equation by parts gives 1
Adding the pressure and velocity equations together and summing over all element D k gives
For the elastic part, taking q = σ , w = v and Theorem 3.
We first consider the case τ u = τ p = 0, which corresponds to a non-dissipative central flux. Then, adding together contributions from integrals on both Ŵ ae and Ŵ ea and consolidating terms involving normal vectors and normal matrices yields
Thus, the contribution from the central portion of the flux sums to zero. Next, we can compute the contribution of penalty
Summing all the contributions, we obtain the desired inequality ∂ ∂t
The stability of the DG formulations (10) and (11)i n Theorem 3.2 requires the use of integration by parts. In order to ensure that this same stability holds at the semi-discrete level, integration by parts must hold when integrals are approximated using quadrature. For affinely mapped simplicial meshes, the geometric terms are constant over each element, such that all spatial integrands on the right-hand side of ( 10) and (11)a r e degree 2N -1 polynomials. Thus, any quadrature which is exact for at least degree 2N -1 polynomials is sufficient for stability.
However, numerous numerical studies demonstrate that, for curved domain boundaries, the use of affinely mapped simplicial meshes limits accuracy to second order [16][17][18][19]. In this section, we assume the triangulation h consists of (possibly curved) elements D k . Under this assumption, the mapping k is no longer affine and the geometric terms are non-constant polynomials within each element. The resulting spatial integrands in (10) and (11)a r e now degree 4N -3 polynomials, while the surface integrands are degree 4N -2 polynomials. Thus, the strength of quadrature required to ensure semi-discrete energy stability of the formulations (10) and (11)i s significantly higher for curved meshes than for affine meshes.
We sidestep these quadrature accuracy requirements on curvilinear meshes by using a "strong-weak" DG formulation, where we discretize the intermediate DG formulation (12)i n Theorem 3.2. Similar formulations have been used to guarantee stability under non-standard basis functions [12,20,21]. Because the formulation ( 12) has already been integrated by parts, the proof of energy stability does not require integrals to be exactly evaluated using quadrature. This quadrature-agnostic stability avoids instability and spurious solution growth for under-integrated DG discretizations on curved meshes [22]. However, it does require an explicit quadrature-based discretization, as opposed to a quadrature-free discretization [19,23].
We outline the matrices involved in a quadrature-based DG discretization in the following section. For simplicity, we now assume constant wavespeed c = 1, such that the strong-weak formulation for the acoustic wave equation is given by
The mass matrix M k is replaced by a weighted mass matrix with weight J k , which we approximate using a weight-adjusted approximation, i.e.
Now, we consider the volume contribution in the pressure equation, i.e.
This contribution becomes more involved to evaluate due to the face that derivatives now lie on the pressure test function q. We follow [2,12] and evaluate this contribution as
where V x q T are quadrature-based differentiation matrices defined by
The terms U x i q are defined at quadrature points as
where x x , ... are evaluations geometric factors at quadrature points and U i denotes the vector of degrees of freedom for the ith velocity component u i . The surface contributions are treated similarly.
In this section, we demonstrate the high order convergence and geometric flexibility of the proposed method. In Section 4.1, we verify that the semi-discrete scheme is energy stable by computing the spectra of the proposed DG schemes. In Section 4.2, we test our method on several classical interface problems with known analytical solutions. In Section 4.3, we implement the proposed scheme on curvilinear meshes and perform convergence analyses. In all numerical experiments, we always choose penalty parameters such that τ u = τ v and τ p = τ σ .
We first verify the energy stability of the proposed method for arbitrary heterogeneous media. We follow the approach in [3] and construct a random stiffness matrix using similarity transforms, such that at every quadrature point, C (x) = UDU T , where D is diagonal matrix with random positive entries d min ≤ D ii ≤ d max and U is a random unitary matrix. For the wavespeed in the acoustic media, we generate positive random values c min ≤ c(x) ≤ c max at quadrature nodes.
Let L denote the matrix induced by the global semi-discrete DG formulation, such that the time evolution of the global solution is governed by with Q denotes a vector of degrees of freedom for (u, p, v, σ ). For practical simulations, taking τ u , τ p > 0r e s u l t s in damping of under-resolved spurious components of the solution.
However, a naive selection of these parameters can result in a more restrictive time-step restriction for stability. We wish to choose τ u τ p as large as possible without increasing the spectra of L when using a central flux (i.e. τ u = τ p = 0). In Fig. 1, we observe that the spectra of L for a central flux is roughly half as large as the spectra of L when taking τ u = τ p = 1. We note that the growth in spectra is due to the large negative real part of the extremal eigenvalues of L, which consistent with the observation that a subset of eigenvalue of L approach -∞ as the penalty parameters increase [24]. Moreover, when we take τ u = τ p = 0.5, the largest real part and imaginary part are most have the same magnitude, which indicates that we can add a dissipative term without shortening the time-step size.
In the following section, we show that the proposed DG method exhibits high order convergence for two classical interface problems: Snell's law and the Scholte wave.
In this experiment, we study the convergence rate of the proposed method for the Snell's law, which models a pressure plane wave incident to an elastic-acoustic interface. The incident wave is reflected as a pressure wave in the acoustic media and transmitted as longitudinal and transverse waves in the elastic media. We follow the problem setting given in [9]. For an incident displacement wave of the form,
The transmitted longitudinal displacement wave is
and the transmitted transverse displacement wave is w ts (x, t) = C ts d ts cos (κ s2 [x 1 sin (α ts ) + x 2 cos (α ts )] -ωt) .
Here, ω is the angular frequency; κ p1 , κ p2 , and κ s2 are wavenumbers of the respective waves and α ip , α rp , α tp and α ts are the associated propagation angles. The displacement directions are
cos (α ts ) sin (α ts ) .
The overall displacement can be written as The wave speeds in each layer are given by
, and the corresponding wavenumbers can be computed from the angular frequency
.
Through Snell's Law, the propagation angles are related to the incident angle α ip
.
The amplitudes of the reflected and transmitted waves are related to the incident wave amplitude
where
cos (α ts ) .
We compute the solution for the specific case of
The computational domain is [-1, 1] 2 and the exact solution is prescribed by tractions on the boundary. Uniform tetrahedral meshes are used in the experiment. Fig. 2 shows the convergence of L 2 errors under mesh refinement for both central fluxes and dissipative penalty fluxes. Optimal O (h N+1 ) rates of convergence are observed for the penalty flux, while an "odd-even" convergence pattern is observed for the central flux.
Scholte waves are boundary waves that propagate along elastic-acoustic interfaces. This problem is designed to the test numerical flux between acoustic and elastic media. In our problem setting, we consider two half-spaces: the upper half, x 2 > 0, is fluid with acoustic material parameters λ 1 , μ 1 = 0, and ρ 1 . The lower half, x 2 < 0, is solid with elastic material parameters λ 2 , μ 2 , and ρ 2 . The displacement of a Scholte wave in the acoustic region is given by
and in the elastic region by
The wavenumber is κ = ω c , with decay rates
, where c is the Scholte wavespeed. The longitudinal and transverse wavespeeds are
.
The wave amplitudes are related to each other through the interface condition (9) 2i 1 -
The Scholte wavespeed c is chosen such that the determinant of ( 14)i s zero, and c satisfies
where r = c/c 2s .
We choose the acoustic and elastic material parameters as λ 1 = 1, ρ 1 = 1, μ 1 = 0, and λ 2 = μ 2 = 1, ρ 2 = 1. For these material parameters, we obtain c = 0.7110017230197 and choose B 1 =-i0.3594499773037, B 2 =-i0.8194642725978, and B 3 = 1. In our experiment, we choose a uniform mesh with different size h covering a square domain [-1, 1] 2 . As with Snell's law, we investigate the convergence rates of the proposed method for a central flux (τ u = τ p = 0) and a penalty flux (τ u = τ p = 1).
Figs. 2 and3 show L 2 error for the Snell's law and Scholte waves at time T = 5, respectively. For penalty fluxes, the computed convergence rate is close to the optimal rate of O (h N+1 ). For central fluxes, we observe again an odd-even pattern, though the rate of convergence is one order lower than observed for Snell's law.
We also computed Scholte wave solutions using more realistic material coefficients from [11]. The fluid media is homogeneous isotropic with an acoustic wavespeed of 1.5 km/s and density 1.0 g/cm 3 . The solid media is homogeneous and isotropic with a P-wave speed of 3.0 km/s and an S-wave speed of 1.5 km/s, with a density of 2.5 g/cm 3 . Errors for a Scholte wave solution at time T = 1a r e shown in Fig. 4.
now numerical experiments verifying the stability and accuracy of the DG scheme presented in Section 3.1 for curvilinear meshes. We use isoparametric mappings in the following experiments, where the mapping from the reference element to each physical element is a polynomial of degree N. We start from a uniform triangular mesh on the square domain =[-1, 1] 2 and place high-order Warp and Blend interpolation nodes on each element. The physical locations (x i , y i ) of these nodes are then perturbed to produce new nodal positions ( x i , y i ), where
x sin (π y) ,
sin (π x) sin (π y) .
These new positions ( x i , y i ) now define a coordinate mapping from the reference element to a curved physical element, producing the warped mesh in Fig. 5(d). This mesh warping is constructed such that x and y deformations of each element are of roughly the same magnitude, while leaving the positions of nodes on the boundary unchanged. both central and penalty fluxes, the real part of all eigenvalues is non-positive (up to machine precision), verifying that the proposed DG scheme is energy stable. The introduction of the curvilinear warping appears to result in a magnification of the real and imaginary parts larger magnitude eigenvalues, which also induces a smaller time-step size.
We compute L 2 errors on a sequence of refined curvilinear meshes for N = 1, 2, 3, 4. From Fig. 6, we observe the rates of convergence of L 2 errors are consistent with the rates observed for affine meshes in Section 4.2.
In this section, we demonstrate the accuracy and flexibility of the proposed DG method for some application-based problems. In the first example, we simulate wave propagation through heterogeneous anisotropic media. In the second example, we present an application of the new DG method to an inverse problem in photoacoustic tomography (PAT).
We examine a model wave propagation problem in heterogeneous anisotropic media. In our experiments, we use two different experimental settings based on [25]. We divide the domain into three parts and set the left half (i.e. x < 0) to be anisotropic elastic media, the right-bottom part (i.e. x > 0, y < 0) to be isotropic elastic media, and the right-upper part (i.e. x > 0, y > 0) to be acoustic media. We assume that the density ρ = 7100 is constant over the whole domain.
In the first experiment, we simulate wave propagation through homogeneous media. The entries of the stiffness matrix C in the anisotropic media are taken to be In the second experiment, we introduce sub-cell heterogeneities to the material parameters. For the isotropic elastic region x < 0, y > 0, we set In the acoustic domain x > 0, y > 0, we again set In all experiments, we set the order of approximation N = 5. We use a uniform triangular mesh of 32768 elements on domain [-0.32, 0.32] 2 . Forcing is applied to the y component of the velocity using a Ricker wavelet point source
where x 0 =-0.02, f 0 = 0.17, and t 0 = 1/ f 0 .
In all implementations, we take the penalty parameters to be τ u = τ p = 1/2. For this value of τ and for the acoustic wave equation in homogeneous media, penalty flux coincides with the upwind flux. Moreover, numerical results suggest that the maximum stable time-step size for τ = 1/2i s the same as the maximum stable time-step size for τ = 0 [ 3], which suggests that this level of dissipation does not require a more restrictive CFL condition. Figs. 7 and8 show the y component of velocity at times T = 30 µs and T 60 µs. In the elastic regions, the results agree with the reference results in [25].
In the elastic-acoustic regions, we observe the presence of propagating pressure wave, while the stress wave ends in a Scholte-type wave propagating along the elastic-acoustic interface. Fig. 8 illustrates the effect of media heterogeneities, which manifest as a spatially-dependent warping of the solution.
Photoacoustic tomography (PAT) is an imaging modality which takes advantage of high-contrast exhibited by optical absorption and the high resolution available for broadband acoustic waves in soft biological tissues. PAT relies on the socalled "photoacoustic effect": a short microwave or light pulse is sent through a patient's body which slightly heats up tissue. The expansion due to heat generates weak acoustic waves, which are measured away from the patient's body. The main step of PAT the recovery of the initial acoustic profile, which in turn provides information about the rate of absorption and tissue properties at different points in the body.
Given the initial state of the pressure field P , the forward mapping F propagates the wave field to the measurements M (Dirichlet data) on the boundary (0, T ) × ∂ . In practice, to produce synthetic measurements, an absorbing boundary condition is employed to allow the waves to radiate outwardly without spurious reflections. The goal of PAT is to invert the forward mapping F : P → M. Typically, a time-reversal method is utilized to approximately invert this forward mapping. The time-reversal mapping R consists of running the wave system backwards in time, from vanishing final condition at {t = T } × , driven from the boundary (0, T ) × ∂ by the time-reversed boundary measurements M as Dirichlet data. The resulting pressure profile at {t = 0} × is an approximation of the original profile P . This approach can be inaccurate for short times and heterogeneous media. However, the quality of the reconstruction can be improved by approximating the exact inversion operator using a truncated Neumann series [26]. Similar reconstruction algorithms have been introduced for several variations of the wave equation [27][28][29][30][31]. We follow the approach proposed in [32], which is summarized in Algorithm 1. These approaches rely on the following error estimate,
which is verifiable when the wave speed is non-trapping (see details in [27]). In other words, the time-reversal mapping R inverts the forward operator F up to a contraction mapping. Algorithm 1 is then the application of a fixed point iteration or truncated Neumann series. The error associated with the nth iteration satisfies,
Algorithm 1 Time-reversal algorithm for PAT. Apply the forward solver with initial condition P n-1 and absorbing boundary conditions. Store the solution at time t = T in P f .
Apply the backward solver with initial condition P f and zero Dirichlet boundary condition. Store the solution at time t = 0i n P b .
Update P n = P n-1 + P b .
We test our PAT algorithm by reconstructing portions of the Shepp-Logan phantom (SLP), which is a standard test for image reconstruction algorithms. The SLP is defined as the sum of 10 ellipses inside the computational domain [-1, 1] 2 .
The specific setting of our experiment is presented in Table 1, and we arbitrarily set the penalty parameters to be τ u = τ p = τ σ = τ v = 1. We simply use the even polynomial function in [33]t o construct a smoothed Shepp-Logan phantom for our numerical simulations with smoothing parameters m = 2, n = 4.
We modify the typical SLP to emulate physical settings found for a human skull. We consider the domain inside domain of Ellipse a and outside of Ellipse b as skull modeled by elastic media. The rest of the domain is acoustic. The meshes (see in Fig. 9 and 12) for the SLP is generated by MESH2D [34], a MATLAB-based mesh-generator for two-dimensional geometries. We use two meshes to test our PAT solver and compare results. The fine mesh consists of 7626 nodes and 14994 elements. The thinnest portion of the elastic domain is resolved using three layers of elements. The coarse mesh consists of 4190 nodes and 8122 elements, and the thinnest portion of the elastic strip is resolved using only one or two layers of elements.
We generate synthetic boundary data by running a forward problem and saving boundary measurements up to final time T = 2. We implement two versions of PAT: the first uses forward and backward solvers based on the discussed elasticacoustic DG formulation, while the second uses a purely acoustic solver for comparison. The wavespeed for the purely acoustic solver is set to be the pressure wavespeed for the elastic system. All experiments are run on an Nvidia TITAN GPU, and the solvers are implemented in the Open Concurrent Compute Abstraction framework (OCCA) [35]f o r clarity and portability.
The relative L 2 errors during each iteration are presented in 2. We observe that, independently of the mesh size, the relative errors of the reconstructed initial data are ≈ 0.06, while the relative errors of the reconstruction from purely acoustic time-reversal are roughly twice as large ≈ 0.12. We present reconstructed initial pressures for both meshes in Fig. 10 and13. From these figures, we observe that using a purely acoustic solver results in larger background noise than using a coupled elastic-acoustic solver. We also observe that the error in the reconstruction is concentrated near the boundary of eclipses and at the elastic-acoustic interfaces (see Figs. 11 and14). The former is due to high gradients in the solution, while the latter may be due to the retention of energy within the elastic region.
In this paper, we present a high order discontinuous Galerkin method for wave propagation in coupled elastic-acoustic media. The method utilizes easily invertible weight-adjusted approximations of weighted mass matrices, as well as an upwind-like penalty numerical flux across the interface between elastic and acoustic media. The formulation is provably discretely energy stable and consistent on arbitrary heterogeneous media, including anisotropy and sub-cell microheterogeneities. An extension of the method to curvilinear meshes achieves similar results. Numerical examples confirm the high order accuracy of this method for analytic solutions to classical interface problems, and results produced by the proposed method are consistent with existing results for isotropic and anisotropic heterogeneous media.
Future work includes the acceleration of the proposed method using tailored Bernstein-Bezier algorithms [36,37], which can reduce the computational complexity of the implementation from O (N 2d ) to O (N d+1 ) in d dimensions, as well as extensions to wave propagation in acoustic-elastic-poroelastic media [38].
The stable and consistent imposition of boundary conditions is described in[1,3].