High Order Finite Difference Methods for the Wave Equation with Non-conforming Grid Interfaces

Abstract

We use high order finite difference methods to solve the wave equation in the second order form. The spatial discretization is performed by finite difference operators satisfying a summation-by-parts property. The focus of this work is on numerical treatments of non-conforming grid interfaces and non-conforming mesh blocks. Interface conditions are imposed weakly by the simultaneous approximation term technique in combination with interface operators, which move discrete solutions between grids at an interface. In particular, we consider an interpolation approach and a projection approach with corresponding operators. A norm-compatible condition of the interface operators leads to energy stability for first order hyperbolic systems. By imposing an additional constraint on the interface operators, we derive an energy estimate of the numerical scheme for the second order wave equation. We carry out eigenvalue analyses to investigate the additional constraint and its relation to stability. In addition, a truncation error analysis is performed, and discussed in relation to convergence properties of the numerical schemes. In the numerical experiments, stability and accuracy properties of the numerical scheme are further explored, and the practical usefulness of non-conforming grid interfaces and mesh blocks is discussed in two practical examples.

Appendices

Appendix 1

We use L, $R^1$ and $R^2$ as subscripts or superscripts to distinguish variables in the left, bottom-right and top-right blocks in the mesh shown in Fig. 2a. The grid functions u, $v_1$ and $v_2$ are the numerical solution vectors in $\varOmega _L$, $\varOmega _{R^1}$ and $\varOmega _{R^2}$, respectively. The mesh blocks are non-conforming, and the grid interfaces can be non-conforming as well. Furthermore, there is no strict restriction on the mesh refinement ratio.

1.1 Semi-discretization

In the following, we construct the penalty terms for the numerical treatment of the T-junction interface conditions. The penalty terms for the outer boundaries are omitted. The semi-discretized equation reads

$$\begin{aligned} \begin{aligned} u_{tt}&=\varvec{D_{2xL}+D_{2yL}}u+SAT_{u1}+SAT_{u2}+SAT_{\partial u}, \\ \begin{pmatrix} v_1 \\ v_2 \end{pmatrix}_{tt}&=\varvec{\begin{pmatrix} D_{2xR^1}+D_{2yR^1} &{}\, \\ \,&{}\quad D_{2xR^2}+D_{2yR^2} \end{pmatrix}}\begin{pmatrix} {v_1} \\ {v_2} \end{pmatrix}+SAT_{v1}+SAT_{v2}+SAT_{\partial v}, \end{aligned} \end{aligned}$$

(20)

where

$$\begin{aligned} \begin{aligned}&SAT_{u1}=\frac{1}{2}\varvec{H_{xL}^{-1}S_{xL}^T}\left( \varvec{E_{0L}}u- \varvec{I_{R2L}}\begin{pmatrix} v_1 \\ v_2 \end{pmatrix} \right) , \\&SAT_{u2}=-\tau \varvec{H_{xL}^{-1}}\left( \varvec{E_{0L}}u- \varvec{I_{R2L}}\begin{pmatrix} v_1 \\ v_2 \end{pmatrix} \right) , \\&SAT_{\partial u}=-\frac{1}{2}\varvec{H_{xL}^{-1}}\left( \varvec{E_{0L}S_{xL}}u-\varvec{I_{R2L}} \begin{pmatrix} \varvec{S_{xR^1}} &{}\, \\ \,&{}\quad \varvec{S_{xR^2}} \end{pmatrix} \begin{pmatrix} v_1 \\ v_2 \end{pmatrix}\right) , \\&\varvec{I_{R2L}}=(E_{0L}\otimes P_{g2f}^L) \varvec{I_{v_1v_2}^u} \begin{pmatrix} E_{LR^1}\otimes P_{f2g}^{R^1} &{}\, \\ \,&{}\quad E_{LR^2}\otimes P_{f2g}^{R^2} \end{pmatrix}, \end{aligned} \end{aligned}$$

(21)

and

$$\begin{aligned} \begin{aligned}&SAT_{v1}=-\frac{1}{2}\begin{pmatrix} \varvec{H_{xR^1}^{-1}S_{xR^1}^T} &{} \\ &{} \varvec{H_{xR^2}^{-1}S_{xR^2}^T} \end{pmatrix}\left( \begin{pmatrix} \varvec{E_{0R^1}} &{}\, \\ \,&{}\quad \varvec{E_{0R^2}} \end{pmatrix} \begin{pmatrix} v_1 \\ v_2 \end{pmatrix} - \varvec{I_{L2R}}u \right) , \\&SAT_{v2}=-\tau \begin{pmatrix} \varvec{H_{xR^1}^{-1}} &{} \\ &{} \varvec{H_{xR^2}^{-1}} \end{pmatrix}\left( \begin{pmatrix} \varvec{E_{0R^1}} &{}\, \\ \,&{}\quad \varvec{E_{0R^2}} \end{pmatrix} \begin{pmatrix} v_1 \\ v_2 \end{pmatrix} - \varvec{I_{L2R}}u \right) , \\&SAT_{\partial v}=\frac{1}{2}\begin{pmatrix} \varvec{H_{xR^1}^{-1}} &{}\, \\ \,&{}\quad \varvec{H_{xR^2}^{-1}} \end{pmatrix}\left( \begin{pmatrix} \varvec{E_{0R^1}S_{xR^1}} &{}\, \\ \,&{}\quad \varvec{E_{0R^2}S_{xR^2}} \end{pmatrix} \begin{pmatrix} v_1 \\ v_2 \end{pmatrix} - \varvec{I_{L2R}S_{xL}}u \right) , \\&\varvec{I_{L2R}}=\begin{pmatrix} E_{R^1L}\otimes P_{g2f}^{R^1} &{}\, \\ \,&{}\quad E_{R^2L}\otimes P_{g2f}^{R^2} \end{pmatrix} \varvec{I_u^{v_1v_2}} (E_{0L}\otimes P_{f2g}^L). \end{aligned} \end{aligned}$$

(22)

The elements in the matrices $\varvec{I_{v_1v_2}^u}$ and $\varvec{I_u^{v_1v_2}}$ are either zero or one. These two matrices are used to pick up the solution along the interface and put it in the right place for coupling.

1.2 Energy Analysis

We note that there is a similarity between the semi-discretization (20) and (3). Comparing the penalty terms in (21) with the penalty terms in (4), the difference is that $E_{LR}\otimes I_{F2C}$ in (4) is replaced by $\varvec{I_{R2L}}$ in (21). Similarly, $E_{RL}\otimes I_{C2F}$ in (5) is replaced by $\varvec{I_{L2R}}$ in (22). As a consequence, the energy analysis for (20) follows in the same way as in Theorem 1, if the corresponding norm-compatible condition (Definition 2) and the norm-contracting condition (Definition 3) are satisfied.

For the T-junction interface in consideration, the norm-compatible condition is

$$\begin{aligned} \begin{pmatrix} H_{yR^1}P_{g2f}^{R^1} &{}\, \\ \,&{}\quad H_{yR^2}P_{g2f}^{R^2} \end{pmatrix}P_{f2g}^L=\left( H_{yL}P_{g2f}^L\begin{pmatrix} P_{f2g}^{R^1} &{} \,\\ \,&{}\quad P_{f2g}^{R^2} \end{pmatrix}\right) ^T, \end{aligned}$$

(23)

and the norm-contracting condition is that the matrices

$$\begin{aligned} \begin{aligned}&H_{yL}\left( I_{yL}-P_{g2f}^L\begin{pmatrix} P_{f2g}^{R^1} &{}\, \\ \,&{}\quad P_{f2g}^{R^2} \end{pmatrix}\begin{pmatrix} P_{g2f}^{R^1} &{}\, \\ \,&{}\quad P_{g2f}^{R^2} \end{pmatrix}P_{f2g}^L\right) ,\\&\begin{pmatrix} H_{yR^1} &{}\, \\ \,&{}\quad H_{yR^2} \end{pmatrix}\left( \begin{pmatrix} I_{yR^1} &{}\, \\ \,&{}\quad I_{yR^2} \end{pmatrix}-\begin{pmatrix} P_{g2f}^{R^1} &{}\, \\ \,&{}\quad P_{g2f}^{R^2} \end{pmatrix}P_{f2g}^L P_{g2f}^L\begin{pmatrix} P_{f2g}^{R^1} &{}\, \\ \,&{}\quad P_{f2g}^{R^2} \end{pmatrix}\right) \end{aligned} \end{aligned}$$

(24)

are symmetric positive semi-definite. It is straightforward to prove (23) by using (15) in Definition 4, so (23) is satisfied for all the projection operators in [14]. However, the validation of the norm-contracting property relies on eigenvalue analyse. We have performed an eigenvalue analysis for the particular mesh setting used in Sect. 5.2.2, and found that the norm-contracting condition is satisfied for the second and fourth order accurate cases. The validation remains with a refined mesh.

Appendix 2

We give a detailed description for the construction of the two grids used in the efficiency studies in Sect. 5.3. A general block $\mathcal {B}_i$ of a decomposition has four boundaries defined by the parametrized curves

$$\begin{aligned} \begin{aligned}&\mathcal {C}_{iS} = \left( x_{iS}(\xi ),y_{iS}(\xi )\right) ,\quad \mathcal {C}_{iN} = \left( x_{iN}(\xi ),y_{iN}(\xi )\right) ,\\&\mathcal {C}_{iW} = \left( x_{iW}(\eta ),y_{iW}(\eta )\right) , \quad \mathcal {C}_{iE} = \left( x_{iE}(\eta ),y_{iE}(\eta )\right) , \end{aligned} \end{aligned}$$

where $0\le \xi \le 1, 0\le \eta \le 1$. $\mathcal {C}_{iS}$ and $\mathcal {C}_{iN}$ describe one pair of opposing sides and $\mathcal {C}_{iW}$ and $\mathcal {C}_{iE}$ the other pair. Let $\mathcal {P}_{iSW}$ denote the point of intersection between the curves $\mathcal {C}_{iS}$ and $\mathcal {C}_{iW}$ e.t.c. A bijection $(x,y) = \mathcal {T}_i(\xi ,\eta )$ from the unit square $\mathcal {S} = [0,1]^2$ to the block $\mathcal {B}_i$ of the decomposition is given by the transfinite interpolation [12] as

$$\begin{aligned} \mathcal {T}_i(\xi ,\eta )= & {} (1-\eta ) \mathcal {C}_{iS} + \eta \mathcal {C}_{iN} + (1-\xi ) \mathcal {C}_{iW} + \xi \mathcal {C}_{iE}\\&-\,\xi \eta \mathcal {P}_{iNE} - \xi (1-\eta ) \mathcal {P}_{iSE} - \eta (1-\xi ) \mathcal {P}_{iNW} - (1-\xi )(1-\eta )\mathcal {P}_{iSW}. \end{aligned}$$

The unit square $\mathcal {S}$ is discretized by the points

$$\begin{aligned} \begin{aligned}&\xi _j = j h_{\xi }, h_{\xi } = 1/(N_{i\xi }-1), \quad j = 0, \dots , N_{i\xi }-1,\\&\eta _k = k h_{\eta }, h_{\eta } = 1/(N_{i\eta }-1), \quad k = 0, \dots , N_{i\eta }-1, \end{aligned} \end{aligned}$$

where $N_{i\xi }$ and $N_{i\eta }$ are integers determining the number of grid points in the spatial directions of the discretization of the block $\mathcal {B}_i$. The corresponding grid points are computed as

$$\begin{aligned} (x_j,y_k) = \mathcal {T}_i(\xi _j,\eta _k). \end{aligned}$$

We now give details on how the grids in $\mathcal {D}_1 - \mathcal {D}_2$ are constructed. $\mathcal {D}_1$ represents a circular cavity in an infinite surrounding medium. We construct the computational grid such that the cavity of radius a is centered inside a square of side 2D. This is done by introducing four blocks $\mathcal {B}^{(1)}_1$–$\mathcal {B}^{(1)}_4$. The bounding curves of block $\mathcal {B}^{(1)}_1$ are given by

$$\begin{aligned} \begin{aligned}&\mathcal {C}^{(1)}_{1S} = a\left( \xi \sqrt{2} - 1/\sqrt{2}, \sqrt{1-(\xi \sqrt{2} - 1/\sqrt{2})^2}\right) ,\\&\mathcal {C}^{(1)}_{1N} = \left( D, 2D \xi - D\right) ,\\&\mathcal {C}^{(1)}_{1W} = \left( -\eta (D-a/\sqrt{2}) - a/\sqrt{2}, \eta (D-a/\sqrt{2}) + a/\sqrt{2}\right) ,\\&\mathcal {C}^{(1)}_{1E} = \left( \eta (D-a/\sqrt{2}) + a/\sqrt{2}, \eta (D-a/\sqrt{2}) + a/\sqrt{2}\right) . \end{aligned} \end{aligned}$$

The bounding curves of the remaining blocks $\mathcal {B}^{(1)}_2$–$\mathcal {B}^{(1)}_4$ that constitute the square with the cavity at the center are obtained via rotation by a factor $\pi /2$,

$$\begin{aligned} \mathcal {C}^{(1)}_{ij} = \mathcal {C}^{(1)}_{i-1j} \begin{bmatrix} \cos {\pi /2}&\quad -\sin {\pi /2}\\ \sin {\pi /2}&\quad \cos {\pi /2}\end{bmatrix},\quad i = 2,3,4,\ \ j = S,N,W,E. \end{aligned}$$

The domain $\mathcal {D}_1$ can now be represented by attaching the grid representing the cavity to one or more Cartesian blocks.

The circular inclusion of $\mathcal {D}_2$ is decomposed into five blocks $\mathcal {B}^{(2)}_1$–$\mathcal {B}^{(2)}_5$. The block $\mathcal {B}^{(2)}_1$ is a square at the center of the circular inclusion with corners at the points $(\pm ad,\pm ad), 0 < d < \sqrt{2}/2$ defined by the bounding curves,

$$\begin{aligned} \begin{aligned}&\mathcal {C}^{(2)}_{1S} = a\left( 2d \xi -d, -d\right) ,\ \mathcal {C}^{(2)}_{1N} = a\left( 2d\xi -d, d\right) ,\\&\mathcal {C}^{(2)}_{1W} = a\left( -d, 2d\eta -d\right) ,\ \mathcal {C}^{(2)}_{1E} = a\left( d, 2d\eta -d\right) . \end{aligned} \end{aligned}$$

Here a is the radius of the circular inclusion. The block $\mathcal {B}^{(2)}_2$ is defined by its bounding curves,

$$\begin{aligned} \begin{aligned}&\mathcal {C}^{(2)}_{2S} = \mathcal {C}^{(2)}_{1N},\\&\mathcal {C}^{(2)}_{2N} = \mathcal {C}^{(1)}_{1S},\\&\mathcal {C}^{(2)}_{2W} = a\left( -\eta (\sqrt{2}/2-d)-d, \eta (\sqrt{2}/2-d) + d\right) ,\\&\mathcal {C}^{(2)}_{2E} = a\left( \eta (\sqrt{2}/2-d) + d, \eta (\sqrt{2}/2-d) + d\right) . \end{aligned} \end{aligned}$$

The bounding curves of the remaining blocks $\mathcal {B}^{(2)}_3$–$\mathcal {B}^{(2)}_5$ that constitute the circular inclusion of $\mathcal {D}_2$ are obtained via rotations by a factor $\pi /2$ as in (6). The inclusion is then centered inside a square of side 2D by attaching it to the blocks $\mathcal {B}^{(1)}_1$–$\mathcal {B}^{(1)}_4$ above. The circular inclusion is now the union of the nine blocks $\mathcal {B}^{(1)}_1$–$\mathcal {B}^{(1)}_4$ and $\mathcal {B}^{(2)}_1$–$\mathcal {B}^{(2)}_5$. The domain $\mathcal {D}_2$ can now be represented by attaching the grid representing the inclusion to one or more Cartesian blocks.