A weight-adjusted discontinuous Galerkin method for wave propagation in coupled elastic-acoustic media

TL;DR

This paper introduces a high-order discontinuous Galerkin method for simulating wave propagation in coupled elastic-acoustic media, emphasizing efficiency, stability, and applicability to complex heterogeneous materials.

Contribution

It develops a novel weight-adjusted DG scheme with an upwind flux for elastic-acoustic interfaces, enhancing efficiency and stability in complex media simulations.

Findings

Method achieves high-order accuracy and energy stability.

Efficient handling of heterogeneous and anisotropic media.

Validated performance in photoacoustic tomography applications.

Abstract

This paper presents a high-order discontinuous Galerkin (DG) scheme for the simulation of wave propagation through coupled elastic-acoustic media. We use a first-order stress-velocity formulation, and derive a simple upwind-like numerical flux which weakly imposes continuity of the normal velocity and traction at elastic-acoustic interfaces. When combined with easily invertible weight-adjusted mass matrices, the resulting method is efficient, consistent, and energy stable on curvilinear meshes and for arbitrary heterogeneous media, including anisotropy and sub-cell (micro) heterogeneities. We numerically verify the high order accuracy and stability of the proposed method, and investigate its performance for applications in photoacoustic tomography.

Tables

Table 1. Table 1: Setting of Shepp-Logan Phantom.

Ellipse	Center	Major Axis	Minor Axis	Theta	Value
a	$(0,0)$	0.69	0.92	0	0
b	$(0,-0.0184)$	0.6624	0.874	0	0
c	$(0.22,0)$	0.11	0.31	$-{0.18}^{\circ }$	0.02
d	$(-0.22,0)$	0.16	0.41	${0.18}^{\circ }$	0.02
e	$(0,0.35)$	0.21	0.25	0	0.01
f	$(0,0.1)$	0.046	0.046	0	0.01
g	$(0,-0.1)$	0.046	0.046	0	0.01
h	$(-0.08,-0.605)$	0.046	0.023	0	0.01
i	$(0,-0.605)$	0.023	0.023	0	0.01
j	$(0.06,-0.605)$	0.023	0.046	0	0.01

Table 2. Table 2: Relative L 2 superscript 𝐿 2 L^{2} errors at each iteration.

Iteration	Fine	Fine (acous)	Coarse	Coarse (acous)
1	0.140530	0.147435	0.140556	0.147103
2	0.094658	0.133881	0.094811	0.133508
3	0.075081	0.130397	0.075347	0.130010
4	0.065585	0.129331	0.065941	0.128939
5	0.060577	0.128973	0.060998	0.128577

Equations

$$\mathbf{x}={\mathbf{\Phi }}^k\hat{\mathbf{x}},\mathbf{x}\in D^k,\hat{\mathbf{x}}\in \hat{D},$$

$$<ltx:note>V_h\left(D^k\right)=V_h\left(\hat{D}\right)\circ {\left({\mathbf{\Phi }}^k\right)}^{-1}=\{{\hat{v}}_h\circ ({\Phi }^k{)}^{-1},{\hat{v}}_h\in V_h(\hat{D}),\},</ltx:note>$$

$$\widehat{D}=\{\left(\widehat{x},\widehat{y},\widehat{z}\right)\geq-1,\quad\widehat{x}+\widehat{y}+\widehat{z}\leq-1\},\qquad V_{h}\left(\widehat{D}\right)=P^{N}\left(\widehat{D}\right)=\big{\{}\widehat{x}^{i}\widehat{y}^{j}\widehat{z}^{k},\quad 0\leq i+j+k\leq N\big{\}}.$$

$$\hat{D}=\{\left(\hat{x},\hat{y}\right)\ge -1,\hat{x}+\hat{y}\le 0\},$$

$$V_{h}\left(\widehat{D}\right)=P^{N}\left(\widehat{D}\right)=\big{\{}\widehat{x}^{i}\widehat{y}^{j},\quad 0\leq i+j\leq N\big{\}}.$$

$$[[p]]=p^{+}-p,[[\mathbf{u}]]={\mathbf{u}}^{+}-\mathbf{u},$$

$$\{\{p\}\}=\frac{1}{2}\left(p^{+}+p\right),\{\{\mathbf{u}\}\}=\frac{1}{2}\left({\mathbf{u}}^{+}+\mathbf{u}\right).$$

$$\begin{array}{r}\frac{1}{\rho c^2}\frac{\partial p}{\partial t}=\nabla \cdot \mathbf{u},\\ \rho \frac{\partial \mathbf{u}}{\partial t}=\nabla p,\end{array}$$

$$0<c_{\text{min}}\le c(\mathbf{x})\le c_{\text{max}}<\infty .$$

$$\begin{array}{rl}{\left(\frac{1}{c^2}\frac{\partial p}{\partial t},q\right)}_{L^2(D^k)}& ={\left(\nabla \cdot \mathbf{u},q\right)}_{L^2(D^k)}+\sum _{f\in \partial D^k}{⟨\frac{1}{2}{\mathbf{n}}^T[[\mathbf{u}]]+\frac{{\tau }_p}{2}[[p]],q⟩}_{L^2(f)},\\ {\left(\frac{\partial \mathbf{u}}{\partial t},\mathbf{w}\right)}_{L^2(D^k)}& ={\left(\nabla p,\mathbf{w}\right)}_{L^2(D^k)}+\sum _{f\in \partial D^k}{⟨\frac{1}{2}[[p]]\mathbf{n}+\frac{{\tau }_u}{2}[[\mathbf{u}]],\mathbf{w}⟩}_{L^2(f)},\end{array}$$

$$\begin{array}{rl}\rho \frac{\partial \mathbf{v}}{\partial t}& =\sum _{i=1}^d{\mathbf{A}}_i^T\frac{\partial \mathbf{\sigma }}{\partial {\mathbf{x}}_i},\\ {\mathbf{C}}^{-1}\frac{\partial \mathbf{\sigma }}{\partial t}& =\sum _{i=1}^d{\mathbf{A}}_i\frac{\partial \mathbf{v}}{\partial {\mathbf{x}}_i},\end{array}$$

$${\mathbf{A}}_1=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right),{\mathbf{A}}_2=\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 1& 0& 0\end{array}\right),$$

$$\begin{array}{rl}{\left(\rho \frac{\partial \mathbf{v}}{\partial t},\mathbf{w}\right)}_{L^2(D^k)}& ={\left(\sum _{i=1}^d{\mathbf{A}}_i^T\frac{\partial \mathbf{\sigma }}{\partial {\mathbf{x}}_i},\mathbf{w}\right)}_{L^2(D^k)}+\sum _{f\in \partial D^k}{⟨\frac{1}{2}{\mathbf{A}}_n^T[[\mathbf{\sigma }]]+\frac{{\tau }_v}{2}{\mathbf{A}}_n^T{\mathbf{A}}_n[[\mathbf{v}]],\mathbf{w}⟩}_{L^2(f)},\\ {\left({\mathbf{C}}^{-1}\frac{\partial \mathbf{\sigma }}{\partial t},\mathbf{q}\right)}_{L^2(D^k)}& ={\left(\sum _{i=1}^d{\mathbf{A}}_i\frac{\partial \mathbf{v}}{\partial {\mathbf{x}}_i},\mathbf{q}\right)}_{L^2(D^k)}+\sum _{f\in \partial D^k}{⟨\frac{1}{2}{\mathbf{A}}_n[[\mathbf{v}]]+\frac{{\tau }_{\sigma }}{2}{\mathbf{A}}_n{\mathbf{A}}_n^T[[\mathbf{\sigma }]],\mathbf{q}⟩}_{L^2(f)},\end{array}$$

$$\begin{array}{rl}{\left({\mathbf{M}}^k\right)}_{ij}& ={\int }_{D^k}{\varphi }_j^k{\varphi }_i^k={\int }_{\hat{D}}{\varphi }_j{\varphi }_i{J}^k,\\ {\left({\mathbf{M}}_w^k\right)}_{ij}& =<ltx:note>{\int }_{D^k}w{\varphi }_j^k{\varphi }_i^k</ltx:note>=<ltx:note>{\int }_{\hat{D}}w{\varphi }_j{\varphi }_i{J}^k</ltx:note>,\\ {\left({\mathbf{M}}_f^k\right)}_{ij}& ={\int }_{\partial D_f^k}{\varphi }_j^k{\varphi }_i^k={\int }_{\partial {\hat{D}}_f}{\varphi }_j{\varphi }_i{J}_f^k,\end{array}$$

$${\left({\mathbf{S}}_1\right)}_{ij}={\int }_{\hat{D}}\frac{\partial {\varphi }_j}{\partial x}{\varphi }_i{J}^k,{\left({\mathbf{S}}_2\right)}_{ij}={\int }_{\hat{D}}\frac{\partial {\varphi }_j}{\partial y}{\varphi }_i{J}^k,{\left({\mathbf{S}}_3\right)}_{ij}={\int }_{\hat{D}}\frac{\partial {\varphi }_j}{\partial z}{\varphi }_i{J}^k.$$

$$\begin{array}{rl}{\mathbf{M}}_{1/c^2}^k\frac{d\mathbf{p}}{dt}& =\sum _{j=1}^d{\mathbf{S}}_j^k{\mathbf{U}}_j+\sum _{f=1}^{N_{\text{faces}}}{\mathbf{M}}_f^k{F}_p\left(\mathbf{p},{\mathbf{p}}^{+},\mathbf{U},{\mathbf{U}}^{+}\right),\\ {\mathbf{M}}^k\frac{d{\mathbf{U}}_i}{dt}& ={\mathbf{S}}_i^k\mathbf{p}+\sum _{f=1}^{N_{\text{faces}}}{\mathbf{n}}_i{\mathbf{M}}_f^k{F}_u\left(\mathbf{p},{\mathbf{p}}^{+},\mathbf{U},{\mathbf{U}}^{+}\right),i=1,\dots ,d,\end{array}$$

$$\begin{array}{rl}{\left({\mathbf{M}}_f^k{F}_p\left(\mathbf{p},{\mathbf{p}}^{+},\mathbf{U},{\mathbf{U}}^{+}\right)\right)}_j& ={\int }_{\partial D_f^k}\frac{1}{2}\left({\tau }_p[[p]]+\mathbf{n}\cdot [[\mathbf{u}]]\right){\varphi }_j^k,\\ {\left({\mathbf{n}}_i{\mathbf{M}}_f^k{F}_u\left(\mathbf{p},{\mathbf{p}}^{+},\mathbf{U},{\mathbf{U}}^{+}\right)\right)}_j& ={\int }_{\partial D_f^k}\frac{1}{2}\left({\tau }_u[[\mathbf{u}]]\cdot \mathbf{n}+[[p]]\right){\varphi }_j^k{\mathbf{n}}_i.\end{array}$$

$$\begin{array}{rl}{\mathbf{M}}_{\rho \mathbf{I}}^k\frac{d\mathbf{V}}{dt}& =\sum _{i=1}^d\left({\mathbf{A}}_i^T\otimes {\mathbf{S}}_i^k\right)\mathbf{\Sigma }+\sum _{f=1}^{N_{\text{faces}}}\left(\mathbf{I}\otimes {\mathbf{M}}_f^k\right)F_v,\\ {\mathbf{M}}_{{\mathbf{C}}^{-1}}^k\frac{d\mathbf{\Sigma }}{dt}& =\sum _{i=1}^d\left({\mathbf{A}}_i\otimes {\mathbf{S}}_i^k\right)\mathbf{V}+\sum _{f=1}^{N_{\text{faces}}}\left(\mathbf{I}\otimes {\mathbf{M}}_f^k\right)F_{\sigma },\end{array}$$

$${\mathbf{M}}_{{\mathbf{C}}^{-1}}^k=\left[\begin{array}{ccc}{\mathbf{M}}_{{\mathbf{C}}_{11}^{-1}}^k& \dots & {\mathbf{M}}_{{\mathbf{C}}_{1d}^{-1}}^k\\ ⋮& \ddots & ⋮\\ {\mathbf{M}}_{{\mathbf{C}}_{d1}^{-1}}^k& \dots & {\mathbf{M}}_{{\mathbf{C}}_{dd}^{-1}}^k\end{array}\right],$$

$$\frac{d\mathbf{U}}{dt}={\left({\mathbf{M}}_{1/c^2}^k\right)}^{-1}{\mathbf{A}}_k\mathbf{U}.$$

$${\mathbf{M}}_w^k\approx {\tilde{\mathbf{M}}}_w^k={\mathbf{M}}^k{\left({\mathbf{M}}_{1/w}^k\right)}^{-1}{\mathbf{M}}^k.$$

$${\left({\mathbf{M}}_w^k\right)}^{-1}\approx {\left({\tilde{\mathbf{M}}}_w^k\right)}^{-1}={\left({\mathbf{M}}^k\right)}^{-1}{\mathbf{M}}_{1/w}^k{\left({\mathbf{M}}^k\right)}^{-1}.$$

$${\mathbf{M}}_{{\mathbf{C}}^{-1}}^{-1}\approx \left(\mathbf{I}\otimes {\mathbf{M}}^{-1}\right){\mathbf{M}}_{\mathbf{C}}\left(\mathbf{I}\otimes {\mathbf{M}}^{-1}\right).$$

$${\left({\mathbf{V}}_q\right)}_{ij}={\varphi }_j\left({\hat{\mathbf{x}}}_i\right),$$

$${\mathbf{M}}^k=J^k\mathbf{M}=J^k{\mathbf{V}}_q^T\text{diag}\left(\hat{\mathbf{w}}\right){\mathbf{V}}_q,{\mathbf{M}}_{c^2}^k=J^k{\mathbf{V}}_q^T\text{diag}\left(\mathbf{d}\right){\mathbf{V}}_q,{\mathbf{d}}_i=\frac{{\hat{\mathbf{w}}}_i}{c^2\left({\mathbf{\Phi }}^k{\hat{\mathbf{x}}}_i\right)}$$

$$\frac{d\mathbf{U}}{dt}={\left({\mathbf{M}}^k\right)}^{-1}{\mathbf{M}}_{c^2}^k{\left({\mathbf{M}}^k\right)}^{-1}{\mathbf{A}}_k\mathbf{U}.$$

$${\left({\mathbf{M}}^k\right)}^{-1}{\mathbf{M}}_{c^2}^k={\mathbf{P}}_q\text{diag}\left(\frac{1}{c^2\left({\mathbf{\Phi }}^k\hat{\mathbf{x}}\right)}\right){\mathbf{V}}_q,$$

$${\mathbf{u}}^{+}\cdot \mathbf{n}=\mathbf{u}\cdot \mathbf{n},p^{+}=p.$$

$${\mathbf{v}}^{+}=\mathbf{v},{\mathbf{A}}_n^T{\mathbf{\sigma }}^{+}={\mathbf{A}}_n^T\mathbf{\sigma }.$$

$$\mathbf{u}\cdot \mathbf{n}=\mathbf{v}\cdot \mathbf{n},{\mathbf{A}}_n^T\mathbf{\sigma }=p\mathbf{n},$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Simulations of wave propagation through elastic-acoustic coupling media are applicable to a ¡ltx:note¿wide¡/ltx:note¿ 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 acoustic-elastic 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 acoustic-elastic media, with elastic-acoustic interface conditions typically incorporated through modifications of the numerical flux. For second order equations, Antonietti et al¡ltx:note¿.¡/ltx:note¿  [7] analyze the stability and convergence of a symmetric interior penalty DG formulation on polygonal and polyhedral meshes. Appelo and Wang [8] introduce 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. ¡ltx:note¿Wilcox et al. [9] derive 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.¡/ltx:note¿ ¡ltx:note¿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.¡/ltx:note¿ Ye et al¡ltx:note¿.¡/ltx:note¿ [11] circumvent 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 ¡ltx:note¿develop¡/ltx:note¿ a high order DG method for acoustic-elastic media based on the first order stress-velocity form of the equations. The proposed method utilizes a simple dissipative upwind-like penalty flux and weight-adjusted 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) ¡ltx:note¿heterogeneities¡/ltx:note¿. ¡ltx:note¿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 acoustic-elastic 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.¡/ltx:note¿

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).

2 Weight-adjusted DG methods for acoustic and elastic wave propagation

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 ¡ltx:note¿Sections 3.1¡/ltx:note¿ and 3.3.

2.1 Mathematical notation

We assume a physical domain $\Omega$, which is exactly represented by a triangulation $\Omega_{h}$ consisting of $K$ non-overlapping elements $D^{k}$. We assume that each element $D^{k}$ is the image of the reference element $\widehat{D}$ under a mapping $\bm{\Phi}^{k}$

[TABLE]

where $\bm{x}=\left(x,y,z\right)$ are physical coordinates on the $k$th element and $\widehat{\bm{x}}=\left(\widehat{x},\widehat{y},\widehat{z}\right)$ are coordinates on the reference element. Over each element $D^{k}$, we define the polynomial approximation space $V_{h}\left(D^{k}\right)$ as

[TABLE]

where $V_{h}\left(\widehat{D}\right)$ is a polynomial approximation space of degree $N$ on the reference element. In this work,111In three dimensions, the reference element and approximation space are the bi-unit right tetrahedron and total degree $N$ polynomials

\widehat{D}=\{\left(\widehat{x},\widehat{y},\widehat{z}\right)\geq-1,\quad\widehat{x}+\widehat{y}+\widehat{z}\leq-1\},\qquad V_{h}\left(\widehat{D}\right)=P^{N}\left(\widehat{D}\right)=\big{\{}\widehat{x}^{i}\widehat{y}^{j}\widehat{z}^{k},\quad 0\leq i+j+k\leq N\big{\}}.

the reference element is taken to be bi-unit right triangle,

[TABLE]

and the reference approximation space $V_{h}\left(\widehat{D}\right)$ is taken to be total degree $N$ polynomials,

[TABLE]

2.2 Discontinuous Galerkin methods for first-order wave equations

On an element $D^{k}$, we define the jump of scalar and vector valued functions across element interfaces as

[TABLE]

where $p^{+},\bm{u}^{+}$ and $p,\bm{u}$ are the neighboring and local traces of the solution over the interface, respectively. ¡ltx:note¿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,+}$.¡/ltx:note¿ The average across an interface is defined as

[TABLE]

In this work, we use a first-order pressure-velocity formulation for the acoustic wave equation (e.g. in fluid media)

[TABLE]

where $p$ is the acoustic pressure, $\bm{u}\in\mathbb{R}^{d}$ is the vector of velocities in each coordinate direction, and $\rho$ and $c$ are density and wavespeed, respectively. For simplicity, we assume unit density $\rho=1$. We also assume that (1) is posed over time $t\in[0,T)$ on the physical domain $\Omega$ with boundary $\partial\Omega$, with the wavespeed bounded from above and below by

[TABLE]

We adopt a DG variational formulation from [12], which is given over element $D^{k}$ by

[TABLE]

where $\bm{n}$ is the ¡ltx:note¿outward¡/ltx:note¿ normal vector, and $\tau_{p},\tau_{u}$ are penalty parameters. Here, $\left(u,v\right)_{L^{2}\left(D^{k}\right)}$ and ${\left\langle u,v\right\rangle}_{L^{2}(f)}$ denote the $L^{2}$ inner products over $D^{k}$ and a face $f$ of the surface $\partial D^{k}$, respectively.

For the elastic wave equation, we use a symmetrized first-order stress-velocity formulation from [3]. Let $\rho$ be the density and $\bm{C}$ be the symmetric matrix form of constitutive tensor relating stress and strain. The first-order system in $d$ dimensions is given by

[TABLE]

where $\bm{v}$ is the vector of velocity and $\bm{\sigma}$ is a vector consisting of unique entries of the symmetric stress tensor. In two dimensions, ¡ltx:note¿the matrices¡/ltx:note¿ $\bm{A}_{i}$ are defined as

[TABLE]

and the expression of $\bm{A}_{i}$ in three dimensions can be found in [3]. ¡ltx:note¿Note that, by factoring out $\bm{C}$, the resulting matrices $A_{i}$ do not involve any material coefficients and all entries are either [math] or $1$.¡/ltx:note¿ The elastic wave equation is discretized using the following DG formulation:

[TABLE]

where $\bm{A}_{n}$ is normal matrix defined as $\bm{A}_{n}=\sum^{d}_{i=1}\bm{n}_{i}\bm{A}_{i}$, and terms $\tau_{v},\tau_{\sigma}$ are penalty parameters introduced on element interfaces. The DG formulations (2) and (4) are provably consistent and energy stable for non-negative penalty parameters $\tau_{p},\tau_{u},\tau_{v},\tau_{\sigma}\geq 0$ [1, 3].

2.3 The semi-discrete matrix system

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}\left(\widehat{D}\right)$ and $V_{h}\left(D^{k}\right)$ are spanned by bases $\{\phi_{i}\}_{i=1}^{N_{p}}$ and $\{\phi_{i}^{k}\}_{i=1}^{N_{p}}$, respectively. The mass matrix $\bm{M}^{k}$, weighted mass matrix $\bm{M}^{k}_{w}$ and face mass matrix $\bm{M}^{k}_{f}$ for $D^{k}$ are defined as

[TABLE]

where $J^{k}$ and $J^{k}_{f}$ are the volume and face Jacobian of the affine mapping $\bm{\Phi}^{k}$, and $w(\bm{x})$ is a spatially varying positive and bounded weight. We also define weak differentiation matrices $\bm{S}_{i}$ with entries

[TABLE]

Using the above notation, the DG formulation (2) can be written in matrix form as

[TABLE]

where $\bm{U}_{i}$ and $\bm{p}$ are degrees of freedom for $\bm{u}_{i}$ and $p$. The flux terms $F_{p},F_{u}$ are defined such that

[TABLE]

The DG scheme (4) for the elastic wave equations can similarly be written as

[TABLE]

where $F_{v},F_{\sigma}$ denote the elastic flux terms, $\otimes$ denotes the Kronecker product, and the matrix-valued weight mass matrix $\bm{M}^{k}_{\bm{C}^{-1}}$ is defined as

[TABLE]

where $\bm{C}^{-1}_{ij}$ denotes the $ij$th entry of $\bm{C}^{-1}$ and $\bm{M}^{k}_{\bm{C}^{-1}_{ij}}$ denotes the scalar weighted mass matrix with weight $\bm{C}^{-1}_{ij}$.

2.4 Weight-adjusted discontinuous Galerkin method

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 $\bm{U}$ denote the vector of all DG degrees of freedom, and let $\bm{M}_{1/c^{2}}^{k},\bm{A}_{k}$ denote the local matrices representing the local DG mass mass matrix and spatial DG formulation, such that the semi-discrete DG scheme can be written over $D^{k}$ as follows:

[TABLE]

When the wavespeed $c^{2}$ is approximated by a constant over each element, it is possible to apply $\left(\bm{M}^{k}_{1/c^{2}}\right)^{-1}$ using only the constant values of $J^{k},c^{2}$ over each element and a single reference mass matrix inverse $\bm{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) is 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 $\widetilde{\bm{M}}^{k}_{w}$

[TABLE]

The inverse of $\widetilde{\bm{M}}^{k}_{w}$ is then

[TABLE]

Since the weight only appears in $\bm{M}^{k}_{1/w}$, $\left(\widetilde{\bm{M}}^{k}_{w}\right)^{-1}$ can be applied using reference inverse mass matrices and a matrix-free quadrature-based evaluation of $\bm{M}^{k}_{1/w}$. Analogously, the inverse of $\bm{M}_{\bm{C}^{-1}}$ can be approximated by the inverse of a matrix-weighted weight-adjusted mass matrix

[TABLE]

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 $\widehat{\bm{x}}_{i},\widehat{\bm{w}}_{i}$ denote the quadrature points and weights on the reference element $\widehat{D}$. We define the interpolation matrix $\bm{V}_{q}$ as

[TABLE]

whose columns consist of values of basis functions at quadrature points. On each element $D^{k}$, we have

[TABLE]

where $\bm{\Phi}^{k}\widehat{\bm{x}}_{i}$ are quadrature points on $D^{k}$ and ¡ltx:note¿ $c^{2}\left(\bm{\Phi}^{k}\widehat{\bm{x}}\right)$ denote the values of the wavespeed at quadrature points¡/ltx:note¿. Plugging the approximation (6) into the local DG formulation (5), we obtain

[TABLE]

Evaluating $\left(\bm{M}^{k}\right)^{-1}\bm{A}_{k}\bm{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) requires applying the product of an unweighted mass matrix and weighted mass matrix. This can be done using quadrature-based matrices as follows:

[TABLE]

where $\bm{P}_{q}=\bm{M}^{-1}\bm{V}_{q}^{T}\textmd{diag}\left(\widehat{\bm{w}}\right)$ is a quadrature discretization of the polynomial $L^{2}$ projection operator on the reference element. Moreover, since $\bm{P}_{q},\bm{V}_{q}$ are reference operators, the implementation of (8) requires only $O\left(N^{d}\right)$ storage for values of the wavespeed $c^{2}\left(\bm{\Phi}^{k}\widehat{\bm{x}}\right)$ at quadrature points for each element. In contrast, storing full weighted mass matrix inverses or factorizations requires $O\left(N^{2d}\right)$ 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 $O(N_{p})\times O(N_{p})$, implying an $O(N^{6})$ storage requirement.

3 Discontinuous Galerkin methods for coupled elastic-acoustic wave equations

For the first-order acoustic and elastic wave equations, the discontinuous Galerkin schemes (2) and (4) are 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 elastic-acoustic coupled 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.

3.1 Upwind-like numerical 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.,

[TABLE]

For an elastic-elastic interface, the velocity and the traction are continuous, i.e.,

[TABLE]

For an interface between elastic and acoustic media, the normal component of the velocity and the traction are continuous, i.e.,

[TABLE]

where $\bm{u}$ and $\bm{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 $\Omega_{e},\Omega_{a}$ denote the elastic and acoustic domains, respectively. Let {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{ea}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{ae}} denote the respective boundaries of $\Omega_{e}$ and $\Omega_{a}$ which correspond to the acoustic-elastic interface. On {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{ae}}, the numerical fluxes are taken to be

[TABLE]

while the numerical fluxes on {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{ea}} are given by

[TABLE]

We now formulate a DG scheme for the first-order elastic-acoustic coupled wave equations. In the acoustic domain $\Omega_{a}$, the DG formulation is given by

[TABLE]

In the elastic domain $\Omega_{e}$, the DG formulation is given by

[TABLE]

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/c^{2},\bm{C}^{-1}$. In our numerical experiments, we approximate the weighted mass matrices induced by micro (sub-cell) heterogeneities in $1/c^{2},\bm{C}^{-1}$ by easily invertible weight-adjusted mass matrices as described in Section 2.4.

3.2 Consistency and energy stability

In this section, we prove that the DG formulations (10) and (11) are consistent and energy stable in arbitrary heterogeneous media.

Theorem 3.1**.**

The coupled discontinuous Galerkin scheme is consistent.

Proof.

Assume that $\bm{u},p,\bm{v},\bm{\sigma}$ are exact solutions of coupled elastic-acoustic wave equations, and that boundary conditions are imposed through consistent modifications of the numerical flux.222The stable and consistent imposition of boundary conditions is described in [1, 3].. Then, plugging them into (10) and (11) causes 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

[TABLE]

At elastic-elastic interfaces, the traction $\bm{A}_{n}^{T}\bm{\sigma}$ and the velocity are continuous, and the numerical flux reduces to

[TABLE]

For an elastic-acoustic interface {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{ae}} ,we have

[TABLE]

Similarly, on {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{ea}}, we have

[TABLE]

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) can also be shown to be energy stable for any choice of $\tau_{u}=\tau_{v}\geq 0,\tau_{p}=\tau_{\sigma}\geq 0$. For simplicity, we assume zero homogeneous Dirichlet boundary conditions on $\partial\Omega$ in the proof of energy stability.

Theorem 3.2**.**

The coupled discontinuous Galerkin scheme is energy stable for $\tau_{u}=\tau_{v}\geq 0,\tau_{p}=\tau_{\sigma}\geq 0$, in the sense that ¡ltx:note¿

[TABLE]

¡/ltx:note¿ where $\Omega_{h}^{a}$ and $\Omega_{h}^{e}$ denote the acoustic and elastic computational domain, respectively.

Proof.

For the acoustic part, taking $q=p,\ \bm{w}=\bm{u}$ and integrating the divergence term of the pressure equation by parts gives

[TABLE]

Adding the pressure and velocity equations together and summing over all element $D^{k}$ gives

[TABLE]

For the elastic part, taking $\bm{q}=\bm{\sigma},\ \bm{w}=\bm{v}$ and Theorem 3.1 in [3] gives

[TABLE]

We first consider the case $\tau_{u}=\tau_{p}=0$, which corresponds to a non-dissipative central flux. Then, adding ¡ltx:note¿together contributions from integrals on both $\Gamma_{ae}$ and $\Gamma_{ea}$¡/ltx:note¿ and consolidating terms involving normal vectors and normal matrices yields

[TABLE]

Thus, the contribution from the central portion of the flux sums to zero. Next, we can compute the contribution of penalty fluxes for $\tau_{u},\tau_{p}>0$

[TABLE]

Summing all the contributions, we obtain the desired inequality

[TABLE]

∎

3.3 Extension to curvilinear meshes

The stability of the DG formulations (10) and (11) in 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) are 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 $\Omega_{h}$ consists of (possibly curved) elements $D^{k}$. Under this assumption, the mapping $\Phi^{k}$ is no longer affine and the geometric terms are non-constant polynomials within each element. The resulting spatial integrands in (10) and (11) are 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) is 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) in 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. ¡ltx:note¿This quadrature-agnostic stability avoids instability and spurious solution growth for under-integrated DG discretizations on curved meshes [22]. ¡/ltx:note¿ 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

[TABLE]

The mass matrix $\bm{M}^{k}$ is replaced by a weighted mass matrix with weight $J^{k}$, which we approximate using a weight-adjusted approximation, i.e.

[TABLE]

Now, we consider the volume contribution in the pressure equation, i.e.

[TABLE]

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

[TABLE]

where $\left(\bm{V}_{q}^{\widehat{x}}\right)^{T}$ are quadrature-based differentiation matrices defined by

[TABLE]

The terms $\bm{U}_{q}^{\widehat{\bm{x}}_{i}}$ are defined at quadrature points as

[TABLE]

where $\bm{x}_{\widehat{x}},\dots$ are evaluations of geometric factors at quadrature points and $\bm{U}_{i}$ denotes the vector of degrees of freedom for the $i$th velocity component $\bm{u}_{i}$. The surface contributions are treated similarly.

4 Numerical experiments

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 $\tau_{u}=\tau_{v}$ and $\tau_{p}=\tau_{\sigma}$.

4.1 Spectra and choice of penalty parameter

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, $\bm{C}\left(\bm{x}\right)=\bm{U}\bm{D}\bm{U}^{T}$, where $\bm{D}$ is diagonal matrix with random positive entries $d_{min}\leq\bm{D}_{ii}\leq d_{max}$ and $\bm{U}$ is a random unitary matrix. For the wavespeed in the acoustic media, we generate positive random values $c_{min}\leq c(\bm{x})\leq c_{max}$ at quadrature nodes.

Let $\bm{L}$ denote the matrix induced by the global semi-discrete DG formulation, such that the time evolution of the global solution is governed by

[TABLE]

with $\bm{Q}$ denotes a vector of degrees of freedom for $(\bm{u},p,\bm{v},\bm{\sigma})$. Figure 1 shows computed eigenvalues of $L$ for different penalty parameters under discretization parameters $N=3$ and $h=1/4$. In both cases, the largest real part of any eigenvalue is $O(10^{-14})$, verifying that the proposed methods are energy stable up to machine precision.

For practical simulations, taking $\tau_{u},\tau_{p}>0$ results 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 $\tau_{u},\tau_{p}$ as large as possible without increasing the spectra of $\bm{L}$ when using a central flux (i.e. $\tau_{u}=\tau_{p}=0$). In Figure 1, we observe that the spectra of $\bm{L}$ for a central flux is roughly half as large as the spectra of $\bm{L}$ when taking $\tau_{u}=\tau_{p}=1$. We note that the growth in spectra is due to the large negative real part of the extremal eigenvalues of $\bm{L}$, which consistent with the observation that a subset of eigenvalue of $\bm{L}$ approach $-\infty$ as the penalty parameters increase[24]. Moreover, when we take $\tau_{u}=\tau_{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.

4.2 Classical interface problems

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.

4.2.1 Snell’s law for an elastic-acoustic interface

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 acoustic-elastic 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,

[TABLE]

the reflected displacement wave is

[TABLE]

The transmitted longitudinal displacement wave is

[TABLE]

and the transmitted transverse displacement wave is

[TABLE]

Here, $\omega$ is the angular frequency; $\kappa_{p1}$, $\kappa_{p2}$, and $\kappa_{s2}$ are wavenumbers of the respective waves and $\alpha_{ip}$, $\alpha_{rp}$, $\alpha_{tp}$ and $\alpha_{ts}$ are the associated propagation angles. The displacement directions are

[TABLE]

The overall displacement can be written as

[TABLE]

The wave speeds in each layer are given by

[TABLE]

and the corresponding wavenumbers can be computed from the angular frequency

[TABLE]

Through Snell’s Law, the propagation angles are related to the incident angle $\alpha_{ip}$

[TABLE]

The amplitudes of the reflected and transmitted waves are related to the incident wave amplitude

[TABLE]

where

[TABLE]

We compute the solution for the specific case of $c_{p1}=1$, $\rho_{1}=1$, $c_{p2}=3$, $c_{s2}=2$, $\rho_{2}=1$, $\omega=2\pi$, $\alpha_{ip}=0.2$, and $C_{ip}=1.0$. 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. Figure 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.

4.2.2 Scholte wave

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 $\lambda_{1}$, $\mu_{1}=0$, and $\rho_{1}$. The lower half, $x_{2}<0$, is solid with elastic material parameters $\lambda_{2}$, $\mu_{2}$, and $\rho_{2}$. The displacement of a Scholte wave in the acoustic region is given by

[TABLE]

and in the elastic region by

[TABLE]

The wavenumber is $\kappa=\frac{\omega}{c}$, with decay rates

[TABLE]

where $c$ is the Scholte wavespeed. The longitudinal and transverse wavespeeds are

[TABLE]

The wave amplitudes are related to each other through the interface condition (9)

[TABLE]

The Scholte wavespeed $c$ is chosen such that the determinant of (14) is zero, and $c$ satisfies

[TABLE]

where $r=c/c_{2s}$.

We choose the acoustic and elastic material parameters as $\lambda_{1}=1,\ \rho_{1}=1,\ \mu_{1}=0$, and $\lambda_{2}=\mu_{2}=1,\ \rho_{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 $\left[-1,1\right]^{2}$. As with Snell’s law, we investigate the convergence rates of the proposed method for a central flux ($\tau_{u}=\tau_{p}=0$) and a penalty flux ($\tau_{u}=\tau_{p}=1$).

Figures 2 and  3 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.

¡ltx:note¿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/cm3. 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/cm3. Errors for a Scholte wave solution at time $T=1$ are shown in Figure 4.¡/ltx:note¿

4.3 Curvilinear meshes

We now present numerical experiments verifying the stability and accuracy of the DG scheme presented in ¡ltx:note¿Section 3.1¡/ltx:note¿ 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 $\Omega=[-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 $(\widetilde{x}_{i},\widetilde{y}_{i})$, where

[TABLE]

These new positions $(\widetilde{x}_{i},\widetilde{y}_{i})$ now define a coordinate mapping from the reference element to a curved physical element, producing the warped mesh in Figure LABEL:fig:curved. 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.

Figure 5 shows computed eigenvalues of the DG discretization matrix for $N=3$ and a warped curvilinear mesh. We use the strong-weak formulation introduced in Section 3.3 and consider both central and penalty fluxes. We observe that for 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 of 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 Figure 6, we observe the rates of convergence of $L^{2}$ errors are consistent with the rates observed for affine meshes in Section 4.2.

5 Application examples

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 and anisotropic media. In the second example, we present an application of the new DG method to an inverse ¡ltx:note¿problem¡/ltx:note¿ in photoacoustic tomography (PAT).

5.1 Heterogeneous anisotropic media

We examine a model wave propagation problem in heterogeneous and 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 $\rho=7100$ is constant over the whole domain.

In the first experiment, we simulate wave propagation through homogeneous media. The entries of the stiffness matrix $\bm{C}$ in the anisotropic media are taken to be

[TABLE]

[TABLE]

and the acoustic wavespeed is set to be

[TABLE]

In the second experiment, we introduce sub-cell ¡ltx:note¿heterogeneities¡/ltx:note¿ to the material parameters. For the isotropic elastic region $x<0,y>0$, we set

[TABLE]

and for the anisotropic elastic region $x<0,y<0$

[TABLE]

In the acoustic domain $x>0,y>0$, we again set

[TABLE]

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

[TABLE]

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 $\tau_{u}=\tau_{p}=1/2$. ¡ltx:note¿For this value of $\tau$ and for the acoustic wave equation in homogeneous media, the penalty flux coincides with the upwind flux. Moreover, numerical results suggest that the maximum stable time-step size for $\tau=1/2$ is the same as the maximum stable time-step size for $\tau=0$ [3], which suggests that this level of dissipation does not require a more restrictive CFL condition¡/ltx:note¿. Figures 7 and 8 show the $y$ component of velocity at times $T=30\mu s$ and $T=60\mu s$. In the elastic regions, the results agree with the reference results in [25]. In the elastic-acoustic regions, we observe the presence of a propagating pressure wave, while the stress wave ends in a Scholte-type wave propagating along the acoustic-elastic interface. Figure 8 illustrates the effect of media heterogeneities, which manifest as a spatially-dependent warping of the solution.

5.2 Photoacoustic tomography

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 so-called “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 is 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 $\mathcal{F}$ propagates the wave field to the measurements $M$ (Dirichlet data) on the boundary $(0,T)\times\partial\Omega$. 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 $\mathcal{F}:P\mapsto M$. Typically, a time-reversal method is utilized to approximately invert this forward mapping. The time-reversal mapping $\mathcal{R}$ consists of running the wave system backwards in time, from vanishing final condition at $\{t=T\}\times\Omega$, driven from the boundary $(0,T)\times\partial\Omega$ by the time-reversed boundary measurements $M$ as Dirichlet data. The resulting pressure profile at $\{t=0\}\times\Omega$ 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,

[TABLE]

which is verifiable when the wave speed is non–trapping (see details in [27]). In other words, the time-reversal mapping $\mathcal{R}$ inverts the forward operator $\mathcal{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 $n^{\text{th}}$ iteration satisfies,

[TABLE]

where $\kappa<1$.

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 $\tau_{u}=\tau_{p}=\tau_{\sigma}=\tau_{v}=1$. We simply use the even polynomial function in [33] to 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 Figure 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 ¡ltx:note¿using¡/ltx:note¿ 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 acoustic-elastic 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] for clarity and portability.

The relative $L^{2}$ errors during each iteration are presented in Table 2. We observe that, independently of the mesh size, the relative errors of the reconstructed initial data are $\approx 0.06$, while the relative errors of the reconstruction from purely acoustic time-reversal are roughly twice as large $\approx 0.12$. We present reconstructed initial pressures for both meshes in Figure 10 and 13. From these figures, we observe that using a purely acoustic solver results in larger background noise than using a coupled acoustic-elastic solver. We also observe that the error in the reconstruction is ¡ltx:note¿concentrated¡/ltx:note¿ near the boundary of eclipses and at the elastic-acoustic interfaces. 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.

6 Conclusion and future work

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 ¡ltx:note¿micro-heterogeneities¡/ltx:note¿. 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].

Acknowledgments

Kaihang Guo and Jesse Chan acknowledge the support of the National Science Foundation under awards DMS-1719818 and DMS-1712639. The work of Sebastian Acosta was partially supported by NSF grant DMS-1712725.