Acoustic Wave Propagation in Complicated Geometries and Heterogeneous Media

Abstract

We construct finite difference discretizations of the acoustic wave equation in complicated geometries and heterogeneous media. Particular emphasis is placed on the accurate treatment of interfaces at which the underlying media parameters have jump discontinuities. Discontinuous media is treated by subdividing the domain into blocks with continuous media. The equation on each block is then discretized with finite difference operators satisfying a summation-by-parts property and patched together via the simultaneous approximation term method. The energy method is used to estimate a semi-norm of the numerical solution in terms of data, showing that the discretization is stable. Numerical experiments in two and three spatial dimensions verifies the accuracy and stability properties of the schemes.

Appendices

Appendix 1: Proof of Lemma 2

The proof of Lemma 2 uses the structure of the $2p$-th order accurate narrow diagonal variable coefficient second derivative SBP operator constructed in [15]. For reference the operators $D_1$ and $D_2^{(b)}$ for the case $p = 2$ are formulated explicitly. The norm $H$ is given by

$$\begin{aligned} H = \varDelta x~\mathrm{diag}\left( \frac{17}{48},\frac{59}{48},\frac{43}{48},\frac{49}{48},1,\dots ,1,\frac{49}{48},\frac{43}{48},\frac{59}{48},\frac{17}{48}\right) . \end{aligned}$$

The 4-th order accurate narrow diagonal first derivative SBP operator $D_1$ is:

$$\begin{aligned} D_1 = \frac{1}{\varDelta x}\begin{bmatrix} \frac{-24}{17}&\quad \frac{59}{34}&\quad \frac{-4}{17}&\quad \frac{-3}{34}&\quad 0&\quad 0\\ \frac{-1}{2}&\quad 0&\quad \frac{1}{2}&\quad 0&\quad 0&\quad 0\\ \frac{4}{43}&\quad \frac{-59}{86}&\quad 0&\quad \frac{59}{86}&\quad \frac{-4}{43}&\quad 0\\ \frac{3}{98}&\quad 0&\quad \frac{-59}{98}&\quad 0&\quad \frac{32}{49}&\quad \frac{-4}{49}\\&\quad&\quad&\quad \frac{1}{12}&\quad \frac{-8}{12}&\quad 0&\quad \frac{8}{12}&\quad \frac{-1}{12}\\&\quad&\quad&\quad&\quad \ddots&\quad \ddots&\quad \ddots&\quad \ddots&\quad \ddots \end{bmatrix} \end{aligned}$$

With the lower right $4 \times 6$ block obtained by rotating the upper left $4 \times 6$ block $180^{\circ }$ and changing the signs of the elements. The part $M^{(b)}$ of the 4-th order accurate narrow diagonal variable coefficient second derivative SBP operator compatible with $D_1$ is

$$\begin{aligned} M^{(b)} = D_1^TH B^{(b)} D_1 + \frac{\varDelta x^5}{18} D_3^TC_3 B_3^{(b)} D_3 + \frac{\varDelta x^7}{144} D_4^TC_4 B_4^{(b)}D_4. \end{aligned}$$

Where

$$\begin{aligned} B^{(b)}&= \mathrm{diag}(b_1,\dots ,b_N) = B^{(b)}_4,\\ B^{(b)}_{3}&= \mathrm{diag}\left( \frac{b_1 + b_{2}}{2},\frac{b_i + b_{i+1}}{2},\dots ,\frac{b_{N-1} + b_{N}}{2}\right) ,\\ C_3&= \mathrm{diag}\left( 0,0,\frac{163928591571}{53268010936},\frac{189284}{185893},1,\dots ,1,\frac{189284}{185893},\frac{163928591571}{53268010936},0,0\right) \\ C_4&= \mathrm{diag}\left( 0,0,\frac{1644330}{301051},\frac{156114}{181507},1,\dots ,1,\frac{156114}{181507},\frac{1644330}{301051},0,0\right) , \\ D_3&= \begin{bmatrix} -1&\quad 3&\quad -3&\quad 1&\quad 0&\quad 0 \\ -1&\quad 3&\quad -3&\quad 1&\quad 0&\quad 0\\ d_{1}&\quad d_{2}&\quad d_{3}&\quad d_{4}&\quad d_{5}&\quad d_{6}\\&\quad&\quad -1&\quad 3&\quad -3&\quad 1\\&\quad&\quad&\quad \ddots&\quad \ddots&\quad \ddots&\quad \ddots \end{bmatrix} \end{aligned}$$

with

$$\begin{aligned} d_{1}&= \frac{-185893}{301051}, d_{2} = \frac{79000249461}{54642863857}, d_{3} = \frac{-33235054191}{54642863857},\\ d_{4}&= \frac{-36887526683}{54642863857}, d_{5} = \frac{26183621850}{54642863857}, d_{6} = \frac{-4386}{181507}. \end{aligned}$$

the lower right $3 \times 6$ block of $D_3$ is obtained by rotating the upper left $3 \times 6$ block $180^{\circ }$ and changing the signs of the elements,

$$\begin{aligned} D_4 = \begin{bmatrix} -1&\quad 3&\quad -3&\quad 1&\quad 0&\quad 0 \\ -1&\quad 3&\quad -3&\quad 1&\quad 0&\quad 0\\ d_{1}&\quad d_{2}&\quad d_{3}&\quad d_{4}&\quad d_{5}&\quad d_{6}\\&\quad&\quad -1&\quad 3&\quad -3&\quad 1\\&\quad&\quad&\quad \ddots&\quad \ddots&\quad \ddots&\quad \ddots \end{bmatrix} \end{aligned}$$

the lower right $3 \times 6$ block of $D_4$ is obtained by rotating the upper left $3 \times 6$ block $180^{\circ }$. The matrix $S$ is given by

$$\begin{aligned} S = \frac{1}{\varDelta _x}\begin{bmatrix} -\frac{11}{16}&\quad 3&\quad -\frac{3}{2}&\quad \frac{1}{3}&\quad 0&\quad 0 \\ 0&\quad 1&\quad 0&\quad 0&\quad 0&\quad 0\\ 0&\quad 0&\quad \ddots&\quad 0&\quad 0&\quad 0\\ 0&\quad 0&\quad 0&\quad 0&\quad 1&\quad 0\\ 0&\quad 0&\quad -\frac{1}{3}&\quad \frac{3}{2}&\quad -3&\quad \frac{11}{16} \end{bmatrix}. \end{aligned}$$

1.1 Proof of Lemma 2

Let $l_p$ be an positive integer, $b_L$ as in the lemma and $B_L$ the matrix with the only non-zero elements $B_{L_{i,i}} = b_L,i = 1\dots l_p$. For $p = 2$ define

$$\begin{aligned} \begin{array}{l} \tilde{B}^{(b)} = B - B_L,\\ \tilde{B}^{(b)}_{3} = B_3 - B_L,\\ \tilde{B}^{(b)}_{4} = B_4 - B_L,\\ \tilde{M}^{(b)} = D_1^T H \tilde{B} D_1 + \frac{\varDelta x^5}{18} D_3^T C_3 \tilde{B}_3 D_3 + \frac{\varDelta x^7}{144} D_4^T C_4 \tilde{B}_4 D_4,\\ M_L = D_1^T H B_L D_1 + \frac{\varDelta x^5}{18} D_3^T C_3 B_L D_3 + \frac{\varDelta x^7}{144} D_4^T C_4 B_L D_4. \end{array} \end{aligned}$$

(32)

By the choice of $B_L$ these matrices are symmetric positive semi-definite. Now, $S^{-T} M_L S^{-1} = \varDelta x b_L \tilde{M}_L$. Where $\tilde{M}_L$ is independent of $\varDelta x$ and $b_L$. By the construction of $B_L$ the only non-zero elements of $\tilde{M}_L$ are located in an upper left square block $A$ of size independent of $\varDelta x$ but dependent on $l_p$, by construction $A = A^T \ge 0$. For $l_p = 1, \dots$ compute $A$ and construct the matrix $B$ with $B_{ij} = A_{ij}, B_{11} = A_{11} - \beta$. $\beta$ is then chosen by computing the eigenvalues of $B$ as the largest number such that $B \ge 0$. Note that also $\beta$ is independent of $\varDelta x$. A choice of $l_p = 4$ gives a value of $\beta$ as in Table 1. Let $\bar{M}_L$ be the matrix resulting from subtracting $\varDelta x b_L \beta$ from the first diagonal element of $\varDelta x b_L \tilde{M}_L$. We get

$$\begin{aligned} \mathbf {u}^T M^{(b)} \mathbf {u} = \mathbf {u}^T (\tilde{M}^{(b)} + M_L) \mathbf {u} = \mathbf {u}^T (\tilde{M}^{(b)} + S^T \bar{M}_L S) \mathbf {u} + \frac{\beta \varDelta x}{b_L} \left( \bar{B}^{(b_L)}S\mathbf {u}\right) ^2_1. \end{aligned}$$

Similarly we can derive the second term of (17) by constructing the corresponding $M_R$ e.t.c., Then taking $\bar{M}^{(b)} = \tilde{M}^{(b)} + M_L + M_R$ proves the lemma for $p = 2$. The proof for $p = 1$ and $p = 3$ follows the same arguments as for $p = 2$. $\square$

Appendix 2: Dirichlet Boundary Conditions

Assume a homogeneous Dirichlet boundary condition at the west boundary. The corresponding SAT term is then

$$\begin{aligned} SAT_{D_{W}} = \varLambda _\rho ^{-1} H^{-1}_x \left( \tau _{D_W} + \sigma _{D} H_y^{-1} \mathbf {F}^{(-)^T}H_y \right) e_{W} \left( \mathbf {u}_{W} - 0 \right) . \end{aligned}$$

(33)

Similar to the discretization of the interface conditions homogeneous Dirichlet boundary conditions adds a term to the total semi-discrete acoustic energy. The following lemma, proved similarly to Lemma 6, determines the parameter $\tau _D$ for this term do be positive semi-definite.

Lemma 7

With the same notation as in Lemma 6 and $\mathbb {1} = \mathrm{diag}(1,\dots ,1)$ the matrix

$$\begin{aligned} \mathcal {T}_D = \begin{bmatrix} -\tau&\quad -\mathbb {1}&\quad -\mathbb {1}\\ -\mathbb {1}&\quad c_1&\quad 0 \\ -\mathbb {1}&\quad 0&\quad c_2 \\ \end{bmatrix} \end{aligned}$$

(34)

is positive semi-definite if $\tau _j \le -\frac{1}{c_{1_j}} - \frac{1}{c_{2_j}}$.

Define

$$\begin{aligned} \mathbf {U}_D = \begin{bmatrix} \mathbf {u}_{W} \\ (B^{(\rho a_{11})} S_x \mathbf {u})_{W} \\ (\varLambda _\rho \varLambda _{\mathbf {a_{12}}} D_{1 y}\mathbf {u})_{W} \end{bmatrix}, H_D = \begin{bmatrix} H&\quad 0&\quad 0\\ 0&\quad H&\quad 0\\ 0&\quad 0&\quad H \end{bmatrix} \end{aligned}$$

(35)

and let

$$\begin{aligned} \tau _{D_j} \le -\frac{1}{b_ {1R_{j,j}}} -\frac{1}{b_ {2R_{j,j}}}, \tau _D = \mathrm{diag}(\tau _{D_1},\dots ,\tau _{D_{N_y}}) \end{aligned}$$

(36)

where $b_ {1R_{j,j}}$ and $b_ {2R_{j,j}}$ are defined by (25), then the addition of the term $\mathbf {U}_D^T H_D \mathcal {T}_D(\tau _D,b_{1 R},b_{2 R}) \mathbf {U}_D$ to the total semi-discrete energy also defines a semi norm of the solution to (21) and (22). Now with a Dirichlet boundary condition. An energy estimate in this semi-norm follows as in Theorem 1.