An improved high order finite difference method for non-conforming grid interfaces for the wave equation

TL;DR

This paper develops a stable, high-order finite difference method for the wave equation on non-conforming grids, extending stability proofs to sixth order accuracy through new penalty terms, and demonstrates improved stability and accuracy.

Contribution

It introduces new penalty terms that enable stability of sixth order schemes on non-conforming interfaces, overcoming previous limitations.

Findings

Sixth order scheme is now provably stable.

Numerical experiments confirm improved accuracy.

Enhanced stability for high-order methods on complex grids.

Abstract

This paper presents an extension of a recently developed high order finite difference method for the wave equation on a grid with non-conforming interfaces. The stability proof of the existing methods relies on the interpolation operators being norm-contracting, which is satisfied by the second and fourth order operators, but not by the sixth order operator. We construct new penalty terms to impose interface conditions such that the stability proof does not require the norm-contracting condition. As a consequence, the sixth order accurate scheme is also provably stable. Numerical experiments demonstrate the improved stability and accuracy property.

Equations

$$u^{T}H{D}_{1}v=-(D_{1}u{)}^{T}Hv-u_0{v}_0+u_N{v}_N,$$

$${\int }_0^{1}U{V}_{x}dx=-{\int }_0^1{U}_{x}Vdx-U(0)V(0)+U(1)V(1),$$

$$M=\tilde{M}+h\theta (BS{)}^{T}BS,$$

$$M^{(b)}={\tilde{M}}^{(b)}+h\sigma b_m(BS{)}^{T}BS,$$

$$b_m=\min (b(x_0),b(x_1),\cdots ,b(x_l),b(x_N),b(x_{N-1}),\cdots ,b(x_{N-l}))$$

$$M^{(b)}=D_1^T{B}^{(b)}H{D}_1+R^{(b)},$$

$$U_{tt}=U_{xx}+U_{yy},(x,y)\in {\Omega }_u,$$

$$V_{tt}=V_{xx}+V_{yy},(x,y)\in {\Omega }_v.$$

$$U(0.5,y,t)=V(0.5,y,t),U_x(0.5,y,t)=V_x(0.5,y,t).$$

$$H_u{I}_{v2u}=(H_v{I}_{u2v}{)}^T.$$

$$H_u(I_u-I_{v2u}{I}_{u2v})\text{ and }{H}_v(I_v-I_{u2v}{I}_{v2u})$$

$$u_{tt}={\mathbf{D}}_{\mathbf{u}}u+SA{T}_{u1}+SA{T}_{u2}+SA{T}_{u3}+SA{T}_{\partial u},$$

$$v_{tt}={\mathbf{D}}_{\mathbf{v}}v+SA{T}_{v1}+SA{T}_{v2}+SA{T}_{v3}+SA{T}_{\partial v},$$

$$SA{T}_{u1}=\frac{1}{2}{\mathbf{H}}_{\mathbf{ux}}^{-1}{\mathbf{S}}_{\mathbf{ux}}^{\mathbf{T}}({\mathbf{E}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u})v),$$

$$SA{T}_{u2}=-\frac{\tau }{2}{\mathbf{H}}_{\mathbf{ux}}^{-1}({\mathbf{E}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u})v),$$

$$SA{T}_{u3}=-\frac{\tau }{2}{\mathbf{H}}_{\mathbf{ux}}^{-1}((E_{ux}\otimes (I_{v2u}{I}_{u2v}))u-(E_{uv}\otimes I_{v2u})v),$$

$$SA{T}_{\partial u}=-\frac{1}{2}{\mathbf{H}}_{\mathbf{ux}}^{-1}({\mathbf{E}}_{\mathbf{ux}}{\mathbf{S}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u}){\mathbf{S}}_{\mathbf{vx}}v),$$

$$SA{T}_{v1}=-\frac{1}{2}{\mathbf{H}}_{\mathbf{vx}}^{-1}{\mathbf{S}}_{\mathbf{vx}}^{\mathbf{T}}({\mathbf{E}}_{\mathbf{vx}}v-(E_{vu}\otimes I_{u2v})u),$$

$$SA{T}_{v2}=-\frac{\tau }{2}{\mathbf{H}}_{\mathbf{vx}}^{-1}({\mathbf{E}}_{\mathbf{vx}}v-(E_{vu}\otimes I_{u2v})u),$$

$$SA{T}_{v3}=-\frac{\tau }{2}{\mathbf{H}}_{\mathbf{vx}}^{-1}((E_{vx}\otimes (I_{u2v}{I}_{v2u}))v-(E_{vu}\otimes I_{u2v})u),$$

$$SA{T}_{\partial v}=\frac{1}{2}{\mathbf{H}}_{\mathbf{vx}}^{-1}({\mathbf{E}}_{\mathbf{vx}}{\mathbf{S}}_{\mathbf{vx}}v-(E_{vu}\otimes I_{u2v}){\mathbf{S}}_{\mathbf{ux}}u).$$

$$\tau \ge \max \left(\frac{1}{2\theta h_{ux}},\frac{1}{2\theta h_{vx}}\right),$$

$$u_t^T(H_{ux}\otimes H_{uy})u_{tt}+v_t^T(H_{vx}\otimes H_{vy})v_{tt}$$

$$\displaystyle=\$$

$$+\frac{1}{2}{u}_t^T{\mathbf{H}}_{\mathbf{uy}}{\mathbf{S}}_{\mathbf{ux}}^{\mathbf{T}}({\mathbf{E}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u})v)$$

$$-\frac{\tau }{2}{u}_t^T{\mathbf{H}}_{\mathbf{uy}}({\mathbf{E}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u})v)$$

$$-\frac{\tau }{2}{u}_t^T{\mathbf{H}}_{\mathbf{uy}}((E_{ux}\otimes (I_{v2u}{I}_{u2v}))u-(E_{uv}\otimes I_{v2u})v)$$

$$-\frac{1}{2}{u}_t^T{\mathbf{H}}_{\mathbf{uy}}({\mathbf{E}}_{\mathbf{ux}}{\mathbf{S}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u}){\mathbf{S}}_{\mathbf{vx}}v)$$

$$-\frac{1}{2}{v}_t^T{\mathbf{H}}_{\mathbf{vy}}{\mathbf{S}}_{\mathbf{vx}}^{\mathbf{T}}({\mathbf{E}}_{\mathbf{vx}}v-(E_{vu}\otimes I_{u2v})u)$$

$$-\frac{\tau }{2}{v}_t^T{\mathbf{H}}_{\mathbf{vy}}({\mathbf{E}}_{\mathbf{vx}}v-(E_{vu}\otimes I_{u2v})u)$$

$$-\frac{\tau }{2}{v}_t^T{\mathbf{H}}_{\mathbf{vy}}((E_{vx}\otimes (I_{u2v}{I}_{v2u}))v-(E_{vu}\otimes I_{u2v})u)$$

$$+\frac{1}{2}{v}_t^T{\mathbf{H}}_{\mathbf{vy}}({\mathbf{E}}_{\mathbf{vx}}{\mathbf{S}}_{\mathbf{vx}}v-(E_{vu}\otimes I_{u2v}){\mathbf{S}}_{\mathbf{ux}}u).$$

$$\frac{d}{dt}\mathbf{G}=\frac{d}{dt}({\mathbf{G}}_1+{\mathbf{G}}_2+{\mathbf{G}}_3)=0,$$

$${\mathbf{G}}_1=$$

$$+u^T(H_{ux}\otimes M_{uy})u+v^T(H_{vx}\otimes M_{vy})v$$

$$+u^T({\tilde{M}}_{ux}\otimes H_{uy})u+v^T({\tilde{M}}_{vx}\otimes H_{vy})v,$$

$${\mathbf{G}}_2=$$

$$-({\mathbf{E}}_{\mathbf{ux}}{\mathbf{S}}_{\mathbf{ux}}u{)}^T{\mathbf{H}}_{\mathbf{uy}}({\mathbf{E}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u})v)$$

$$+\frac{\tau }{2}({\mathbf{E}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u})v{)}^T{\mathbf{H}}_{\mathbf{uy}}({\mathbf{E}}_{\mathbf{ux}}u-(E_{uv}\otimes I_{v2u})v),$$

$${\mathbf{G}}_3=$$

$$-({\mathbf{E}}_{\mathbf{vx}}{\mathbf{S}}_{\mathbf{vx}}v{)}^T{\mathbf{H}}_{\mathbf{vy}}((E_{vu}\otimes I_{u2v})u-{\mathbf{E}}_{\mathbf{vx}}v)$$

$$+\frac{\tau }{2}((E_{vu}\otimes I_{u2v})u-{\mathbf{E}}_{\mathbf{vx}}v{)}^T{\mathbf{H}}_{\mathbf{vy}}((E_{vu}\otimes I_{u2v})u-{\mathbf{E}}_{\mathbf{vx}}v).$$

$$\begin{array}{rl}{J}_U{U}_{tt}& =(a{U}_x{)}_x+(c{U}_y{)}_y+(b{U}_y{)}_x+(b{U}_x{)}_y,\\ J_V{V}_{tt}& =(\alpha V_x{)}_x+(\gamma V_y{)}_y+(\beta V_y{)}_x+(\beta V_x{)}_y,\end{array}$$

Full text

1 Introduction

Wave propagation can be modeled by hyperbolic partial differential equations (PDEs). When solving a hyperbolic PDE by a finite difference method, to achieve a certain accuracy a minimum number of grid points per wavelength is required. This number is smaller with a high order method than with a low order method, which makes high order finite difference methods more efficient to solve wave propagation problems on smooth domains, see the pioneering paper [17] for first order hyperbolic PDEs, and the recent work [13] for second order hyperbolic PDEs.

On a uniform grid, high order central finite difference stencils are easily constructed by Taylor expansions [9]. Close to boundaries, one–sided stencils can be used. It is important that the boundary closure is accurate, and the numerical scheme is stable. One successful approach is to use a finite difference operator satisfying a summation–by–parts (SBP) property [19]. When an SBP operator acts on a grid function, it mimics the integration–by–parts principle in the continuous setting. An energy estimate can be obtained if boundary conditions are imposed appropriately, for example by using the simultaneous–approximation–term (SAT) method [4] or ghost points [29]. A scheme satisfying an energy estimate is called energy stable [11, 12].

Finite difference methods in the basic form can only be used on box–shaped domains. When complex geometry is present, the domain can be partitioned into blocks to resolve the geometrical feature. Each block has four sides and is mapped to a reference domain. If the corners of adjacent blocks meet, we say they are conforming blocks; otherwise they are non–conforming. In addition, a grid interface is conforming if no hanging nodes are present. When partitioning a domain, we can always make the blocks and interfaces conforming. However, in many situations it is desirable to use a more flexible strategy of partition that leads to non–conforming blocks and grid interfaces.

As an example, we consider a wave traveling in a heterogeneous medium with the wave speed varying in space. The wavelength is proportional to the wave speed for a given frequency. For accuracy the grid spacing is determined by the shortest wavelength. If a uniform grid is used in the entire computational domain, then the grid spacing must be small enough to resolve the shortest wavelength, resulting in an unnecessarily fine grid elsewhere. It is then more efficient to construct a grid according to the wavelength in each block, which leads to non–conforming interfaces with hanging nodes.

If only conforming blocks are used, the domain partitioning may end up with many blocks of small size. To use a high order finite difference method, a minimum number of grid points is required in each block due to the stencil width. This then results in unnecessarily many grid points in the small blocks, and consequently a suboptimal performance of the numerical scheme. In such a situation, non–conforming blocks are more appropriate.

In an SBP finite difference method, interface conditions can also be imposed by the SAT method [5, 6] or ghost points [29]. In the SBP–SAT framework, wave propagation in a heterogeneous medium with complex geometry is considered in [33]. A stable and accurate multi–block finite difference method with conforming grid interfaces and blocks is presented. The focus in [36] is the numerical treatment of non–conforming interfaces and blocks by using SBP–preserving interpolation operators [16, 22]. Energy stability is proved with an assumption that the interpolation operators are norm–contracting. In the same paper [36] , it is verified that not all interpolation operators satisfy this assumption, and instability occurs when the sixth order method is used on a domain with non–conforming, curved interfaces.

In this paper, we construct new penalty terms in the SBP–SAT finite difference framework for the numerical interface treatment. The resulting scheme is energy stable even when the interpolation operators are not norm–contracting. This extends the provably stable scheme from fourth order accuracy [36] to sixth order accuracy. The technique can be potentially used to construct even higher order schemes, provided that the corresponding SBP operators exist. Another contribution of this paper is the numerical treatment of non–conforming blocks and interfaces on curvilinear grids, where as in [36] such a case is studied on Cartesian grids. We also conduct numerical experiments to verify that the new sixth order scheme is stable with non–conforming, curved interfaces.

The paper is organized as follows. In Section 2, we introduce the SBP–SAT finite difference method. In Section 3, we consider the wave equation on a Cartesian grid and present the new penalty terms for numerical interface treatments. Stability is proved by the energy method. We then generalize the scheme to non–conforming blocks and grid interfaces on curvilinear grids in Section 4. Numerical experiments are presented in Section 5 to verify the stability and accuracy property of the developed scheme. We draw conclusion in Section 6.

2 SBP–SAT finite difference methods

Finite difference operators satisfying an SBP property have been widely used to discretize time dependent PDEs. An SBP operator has central finite difference stencils in the interior, and special one–sided stencils at a few grid points near boundaries. The boundary stencils are chosen so that the operator satisfies a summation–by–parts property, which is the discrete counterpart of the integration–by–parts principle. With the SAT method imposing boundary and interface conditions, the SBP–SAT finite difference method possesses a great advantage: it is possible to prove energy stability for high order accurate schemes for initial–boundary–value problems.

To introduce the SBP–SAT finite difference method, we consider the one dimensional domain $[0,1]$ discretized by the grid points $x_{j}=jh,j=0,1,\cdots,N$ with a constant grid spacing $h=1/N$. We use the capital letter, for example, $U$, to denote a smooth function in $[0,1]$, and the corresponding small letter, $u$, to denote its values on the grid $u=[U(x_{0}),U(x_{1}),\cdots,U(x_{N})]^{T}$.

2.1 Definitions of SBP operators

The SBP concept and the first derivative SBP operator $D_{1}\approx\partial/\partial x$ are introduced in [19], and later refined in [31]. Formally it is defined as follows.

Definition 1**.**

A difference operator $D_{1}=H^{-1}Q$ approximating $\partial/\partial x$ is a diagonal norm first derivative SBP operator if $H$ is diagonal positive definite and $Q+Q^{T}=\text{diag}(-1,0,\cdots,0,1)$.

The operator $H$ defines the SBP norm, and leads to the identity

[TABLE]

which is the discrete analogue of the integration–by–parts formula

[TABLE]

since the norm $H$ is also a quadrature [7, 14].

For the second derivative, we distinguish between a constant coefficient operator $D_{2}\approx\partial^{2}/\partial x^{2}$ and a variable coefficient operator $D_{2}^{(b)}\approx\partial/\partial x(b(x)\partial/\partial x)$ with a known function $b(x)>0$.

Definition 2**.**

A difference operator $D_{2}=H^{-1}(-M+BS)$ approximating $\partial^{2}/\partial x^{2}$ is a diagonal norm second derivative SBP operator if $H$ is diagonal positive definite, $M$ is symmetric positive semi–definite, $B=\text{diag}(-1,0,\cdots,0,1)$, and the first and last row of $S$ approximate $\partial/\partial x$ at the two boundaries, respectively.

Such an operator is constructed in [25]. It is later found in [3, 23] that the operator $M$ in $D_{2}$ satisfies the following property.

Lemma 1**.**

The symmetric positive semi–definite operator $M$ can be written as

[TABLE]

where $\widetilde{M}$ is also symmetric positive semi–definite, $\theta>0$ is a constant independent of $h$, $B$ and $S$ are the same as in Definition 2.

Lemma 1 is often referred to as the borrowing trick, as we can borrow from the positive semi–definite operator $M$ a small, mesh dependent amount, with the resulting operator $\widetilde{M}$ still positive semi–definite. This property is essential for energy stability of problems with interfaces or Dirichlet boundary conditions.

For the variable coefficient case we have correspondingly

Definition 3**.**

A difference operator $D_{2}^{(b)}=H^{-1}(-M^{(b)}+B^{(b)}S)$ approximating $\partial/\partial x(b(x)\partial/\partial x)$ is a diagonal norm second derivative variable coefficient SBP operator if $H$ is diagonal positive definite, $M^{(b)}$ is symmetric positive semi–definite, $B^{(b)}=\text{diag}(-b(x_{0}),0,\cdots,0,b(x_{N}))$, and the first and last row of $S$ approximate $\partial/\partial x$ at the two boundaries, respectively.

Such an operator is constructed in [22], and the operator $M^{(b)}$ has the following two important properties [33].

Lemma 2**.**

The symmetric positive semi–definite operator $M^{(b)}$ can be written as

[TABLE]

where $\widetilde{M}$ is also symmetric positive semi–definite, $\sigma>0$ is a constant independent of $h$, $B$ and $S$ are the same as in Definition 2, and

[TABLE]

with a constant $l$ independent of $h$.

Lemma 2 for the variable coefficient SBP operators is an analogue of Lemma 1 for the constant coefficient case. We note that $b_{m}$ is the smallest value of the variable coefficient $b(x)$ on the first and last $l$ grid points. The smaller $b_{m}$ is, the less we can borrow from $M^{(b)}$.

Lemma 3**.**

The SBP operator $D_{2}^{(b)}$ is compatible with $D_{1}$ if $M^{(b)}$ can be written as

[TABLE]

where $R^{(b)}$ is symmetric positive semi–definite, and $B^{(b)}$ is the same as in Definition 3.

Lemma 3 is essential for energy stability when mixed derivatives are present in the equation, for example the wave equation on curvilinear grids and the elastic wave equation.

The definitions and precise forms of the above operators can be found in [19, 21, 22, 25, 31]. These operators have the minimal interior stencil width. In addition, they have the same associated norm $H$ for a given accuracy order. In the stability analysis, we only consider numerical treatment of interface conditions. As a consequence, the operator $B$ in the preceding lemmas only has one nonzero element, corresponding to the terms on the interface.

The interior stencil of an SBP operator is the standard central finite difference stencil with truncation error $\mathcal{O}(h^{2p})$. On a few grid points near boundaries, special one-sided stencils are used to fulfill the SBP requirement with a larger truncation error $\mathcal{O}(h^{p})$. Operators $D_{1}$ and $D_{2}$ with $p=1,2,3,4$ are constructed in [19, 31] and [25], respectively. The variable coefficient operators $D_{2}^{(b)}$ with $p=1,2,3$ are constructed in [21].

In this paper, we call the above SBP operators $2p^{th}$* order accurate*. When using in a numerical scheme, we also call the scheme $2p^{th}$* order accurate*, even though the truncation error of the numerical scheme may not be $\mathcal{O}(h^{2p})$ or $\mathcal{O}(h^{p})$. In the discussion of accuracy, we make the truncation error of the scheme precise.

2.2 The SAT method

An SBP operator only approximates a certain derivative, but does not impose any boundary condition. The boundary conditions must be imposed carefully so that an energy estimate can be obtained to ensure stability. This can be done by for example the SAT method [4], the projection method [26, 27, 28] and the ghost points method [29]. In this paper, we choose the SAT method to impose both boundary and interface conditions, since in many cases it is easy to derive an energy estimate. The key ingredient of the SAT method is to add penalty terms to the semi–discretized equation and choose penalty parameters so that an energy estimate is obtained. This technique bears a similarity to the Nitsche’s finite element method [30], and the discontinuous Galerkin method [2, 10]. Detailed discussions of the SBP–SAT finite difference methods can be found in [8, 32].

3 The wave equation on a Cartesian grid

We start by considering the wave equation in two space dimensions in a composite domain $\Omega=[0,1]^{2}$ with an interface $\Gamma$ at $x=0.5$. The left and right domain are denoted by $\Omega_{u}$ and $\Omega_{v}$, respectively, and the equations are

[TABLE]

At the interface the physical conditions are

[TABLE]

For a wellposed problem, suitable boundary conditions must be imposed at the boundaries. As the focus in this paper is the numerical treatment of interface coupling, we exclude discussions on boundary conditions and the corresponding numerical techniques. We refer to [18] for physical boundary conditions, and [24] for the numerical techniques.

To solve (2)–(4), we start by generating a Cartesian grid in each domain independently with $n_{ux}\times n_{uy}$ grid points in $\Omega_{u}$ and $n_{vx}\times n_{vy}$ grid points in $\Omega_{v}$. We are particularly interested in a non–conforming interface when $n_{uy}\neq n_{vy}$. In this case, the solutions on the interface must be interpolated. We denote $I_{u2v}$ and $I_{v2u}$ interpolation operators that interpolate the solution from $\Omega_{u}$ to $\Omega_{v}$, and from $\Omega_{v}$ to $\Omega_{u}$, respectively. In the SBP–SAT finite difference framework, these operators must satisfy certain conditions so that the scheme could be energy stable.

Definition 4**.**

Let $H_{u}$ and $H_{v}$ denote the SBP norms on the interface for the grid in $\Omega_{u}$ and $\Omega_{v}$, respectively. The interpolation operators $I_{u2v}$ and $I_{v2u}$ are norm–compatible if

[TABLE]

In [36], it is also defined that the interpolation operators are norm–contracting if the two operators

[TABLE]

are symmetric positive semi–definite, where $I_{u}$ and $I_{v}$ are identity operators.

Norm–compatible interpolation operators are first constructed in [22] for the case of a 1:2 mesh refinement ratio, and are extended to an arbitrary ratio in [16]. The accuracy property of these interpolation operators has a similar fashion as the corresponding SBP operators. More precisely, the interpolation error is $\mathcal{O}(h^{2p})$ in the interior of the interface, and $\mathcal{O}(h^{p})$ on a few grid points near the edge of the interface. Therefore, the interpolation is exact only for polynomials of order up to $p-1$. In [20], it is proved that it is not possible to construct norm–compatible interpolation operators $I_{u2v}$ and $I_{v2u}$ such that both interpolate polynomials of order $p$ or higher.

A stable SBP–SAT finite difference method for solving (2)–(4) is presented in [36]. Energy stability is proved by assuming the interpolation operators are norm–compatible and norm–contracting. While the norm–compatible condition can be constrained when constructing the operators, it is not easy to take into account the norm–contracting condition. In fact, the interpolation operators with higher than fourth order accuracy in [16, 22] are not norm–contracting.

Below we present a new way of imposing the interface conditions (4) with the advantage that an energy estimate is obtained without requiring the interpolation operators to be norm–contracting. For cleaner notations, terms imposing boundary conditions are omitted.

Equation (2)–(4) are discretized in space as

[TABLE]

where

[TABLE]

and

[TABLE]

The numerical solution vectors $u$ and $v$ approximate the true solution $U$ and $V$, respectively. The solution vectors are arranged column–wise, i.e. the first few elements of $u$ and $v$ correspond to the solutions on the left boundary of $\Omega_{u}$ and $\Omega_{v}$, respectively. Most operators in two space dimensions can be extended from the corresponding one dimensional operators by using a Kronecker product $\otimes$. Such two dimensional operators are denoted by bold letters, with the subscript indicating the spatial direction and the grid function that the operator is associated to. For example, the operator $\boldsymbol{H_{ux}^{-1}}$ equals to $H_{ux}^{-1}\otimes I_{uy}$, where $H_{ux}^{-1}$ is the inverse of the SBP norm in the $x$–direction acting on $u$, and $I_{uy}$ is an identity operator. The operator $\boldsymbol{E}$ extracts the numerical solution at the interface. $\boldsymbol{D_{u}}u$ and $\boldsymbol{D_{v}}v$ are SBP approximations of $U_{xx}+U_{yy}$ and $V_{xx}+V_{yy}$, respectively.

We compare the above scheme with the ones described in [33, 36] by discussing the penalty terms (8a)–(8d). A term like $\boldsymbol{E_{ux}}u$ used in [33, 36] is broken into two parts: $1/2\boldsymbol{E_{ux}}u$ in (8b) and $1/2E_{ux}\otimes I_{v2u}I_{u2v}u$ in (8c). Note the relation between them: $E_{ux}\otimes I_{v2u}I_{u2v}u$ is just $\boldsymbol{E_{ux}}u$ interpolated to the grid on the interface of $\Omega_{v}$, then interpolated back to the grid of $\Omega_{u}$. Since the interpolation is not exact, $E_{ux}\otimes I_{v2u}I_{u2v}u$ differs from $\boldsymbol{E_{ux}}u$ by the truncation error of the interpolation operators. It is this change in penalty terms that makes the scheme stable without requiring the interpolation operators to be norm–contracting. We summarize the stability result in the following theorem.

Theorem 1**.**

With norm–compatible interpolation operators, the semi–discretization (6)–(7) is stable for any $\tau$ such that

[TABLE]

where $\theta$ is the constant in Lemma 1, and $h_{ux}$ and $h_{vx}$ are the mesh size in the $x$–direction in $\Omega_{u}$ and $\Omega_{v}$, respectively.

Proof.

We prove stability by the energy method. Multiplying from the left of (6) by $u_{t}^{T}(H_{ux}\otimes H_{uy})$ and (7) by $v_{t}^{T}(H_{vx}\otimes H_{vy})$, we obtain

[TABLE]

We note that the left–hand side of the above equation can be written as the time derivative of a quadratic term. The main idea of deriving an energy estimate is to move all terms on the right–hand side to the left, and determine the penalty parameter $\tau$ so that all terms on the left–hand side can be written as the time derivative of a non–negative quantity, i.e. the discrete energy.

To do so, we use the borrowing trick in Lemma 1 for the SBP operators in the $x$–direction, and the norm–compatible property (5) of the interpolation operators, to obtain the change of energy $\mathbf{G}$ as

[TABLE]

where

[TABLE]

Clearly, $\mathbf{G_{1}}\geq 0$. By Young’s inequality, we have $\mathbf{G_{2}}\geq 0$ and $\mathbf{G_{3}}\geq 0$ if $\tau\geq 1/(2\theta h_{ux})$ and $\tau\geq 1/(2\theta h_{vx})$, respectively. Therefore, the energy is conserved and the scheme is stable when (10) is satisfied. ∎

We note that in the scheme developed in [36] the energy is greater or equal to $\mathbf{G}$ in (11), with the inequality resulted from the norm–contracting condition.

4 The wave equation on curvilinear grids

In this section, we generalize the scheme to problems on curvilinear grids. We consider two cases: conforming blocks and non–conforming blocks, which are illustrated in Figure 1a and 1b, respectively.

4.1 Numerical interface treatment of conforming blocks

With only conforming blocks in the domain, the corners of adjacent blocks meet. We consider again the domain $\Omega=[0,1]^{2}$ but partitioned into two blocks $\Omega_{u}$ and $\Omega_{v}$ by a curved interface. The grids are then constructed independently in each block, see an illustration in Figure 1a. The grids in each block are mapped to a Cartesian grid in a reference domain. The governing equations are also transformed from the physical domain to the reference domain, and the computation is performed in the reference domain. We refer to the textbook [15] for a detailed discussion on grid generation.

The transformed equation in the reference domain can be derived by using the chain rule in calculus. We omit its derivation, and refer to [1]. With $(x,y)$ denoting the coordinate in the reference domain, the unknown variables $U$ and $V$ are governed by the equation

[TABLE]

where the Jacobians $J_{U},J_{V}>0$ and the variable coefficients satisfy $ac-b^{2}>0$ and $\alpha\gamma-\beta^{2}>0$. The interface conditions become

[TABLE]

on $(x,y)\in\Gamma$. The coefficients in (12) consist of metric derivatives that depend on the geometry of the physical domain and the transformation. The metric derivatives can either be computed analytically, or approximated to sufficient high accuracy. In the experiments in this paper, we choose the latter approach by using a tenth order finite difference stencil, making sure that the approximation of metric terms does not affect the overall accuracy of the numerical scheme.

As shown in [33], the equation (12) together with the interface condition (13) admit a continuous energy estimate, thanks to the positive definiteness of the matrices $\begin{bmatrix}a&b\\ b&c\end{bmatrix}$ and $\begin{bmatrix}\alpha&\beta\\ \beta&\gamma\end{bmatrix}$ because of $ac-b^{2}>0$ and $\alpha\gamma-\beta^{2}>0$.

The equations in (12) are discretized by using the SBP operators, and the two blocks are patched together by the SAT method. The semi–discretized equations are

[TABLE]

where

[TABLE]

and

[TABLE]

We now clarify the notations in the semi–discretization (14).

The mixed–derivative terms: $(bU_{y})_{x}$ is approximated by $\boldsymbol{D_{1ux}}{\Lambda_{b}}\boldsymbol{D_{1uy}}u$, where ${\Lambda_{b}}$ is a diagonal matrix with diagonal entries $b(x,y)$ evaluated on the grid. The operators approximating the other three mixed–derivative terms are constructed in a similar way. 2. 2.

The variable–coefficient terms: In general the variable coefficients are functions of both $x$ and $y$, therefore an operator approximating $(aU_{x})_{x}$ cannot be constructed by a single Kronecker product, but as a sum

[TABLE]

where $a_{i}$ is $a(x_{i},y)$ evaluated on the grid, and $E_{uy}^{i}$ has value one in entry $(i,i)$ and zeros elsewhere. The second derivative operator $D_{2ux}^{(a_{i})}$ is defined in Definition 3 with the operator $M^{(a_{i})}$ satisfying Lemma 2 and 3. The operators approximating the other three variable coefficient terms are constructed in a similar way. 3. 3.

The penalty terms: The two interpolation operators $I_{u2v}$ and $I_{v2u}$ constructed in [16, 22] satisfy the norm–compatible condition

[TABLE]

Similar to the continuous case, the matrices $\Lambda_{u}=\begin{bmatrix}\Lambda_{a}&\Lambda_{b}\\ \Lambda_{b}&\Lambda_{c}\end{bmatrix}$ and $\Lambda_{v}=\begin{bmatrix}\Lambda_{\alpha}&\Lambda_{\beta}\\ \Lambda_{\beta}&\Lambda_{\gamma}\end{bmatrix}$ are positive definite. In fact, the eigenvalues of $\Lambda_{u}$ and $\Lambda_{v}$ play an important role in the stability analysis. In particular, we will use the smallest eigenvalue

[TABLE]

where $i=1,2,\cdots,n_{ux},\ j=1,2,\cdots,n_{uy},\ k=1,2,\cdots,n_{vx},\ l=1,2,\cdots,n_{vy}$. Note that $\delta>0$, $a_{ij}-\delta\geq 0$ and $\alpha_{ij}-\delta\geq 0$.

Let $a_{\max}$, $b_{\max}$, $\alpha_{\max}$ and $\beta_{\max}$ denote the maximum values of the variable coefficients $a(x,y)$, $b(x,y)$, $\alpha(x,y)$ and $\beta(x,y)$ evaluated on the interface, respectively. We also denote $h_{ux}$ and $h_{vx}$ the mesh size in the $x$–direction in $\Omega_{u}$ and $\Omega_{v}$, respectively. Stability of the semi–discretization is given by the following theorem.

Theorem 2**.**

The semi–discretization (14)–(16) is stable if the interpolation operators $I_{u2v}$ and $I_{v2u}$ are norm–compatible and the penalty parameter $\tau$ satisfies

[TABLE]

where $\sigma$ is defined in Lemma 2, and $\delta$ is defined in (17).

Proof.

See Appendix. ∎

We remark that it is possible to use a penalty parameter that varies on grid points, and such a scheme is proposed in [33] for problems with conforming interfaces. The penalty parameter $\tau$ used in Theorem 2 corresponds to the largest value on all grid points. We also note that when the penalty parameter is chosen to be equal to the stability limit, accuracy deduction has been observed, and proved in some settings by a normal mode analysis [34]. Though there is no upper bound of $\tau$ for energy stability, a very large penalty parameter increases the spectral radius of the spatial discretization, and leads to a small time step.

For accuracy, the SBP operators have truncation error $\mathcal{O}(h^{2p})$ in the interior and $\mathcal{O}(h^{p})$ near the boundaries and the interface. The interpolation operators constructed in [16, 22] have truncation error $\mathcal{O}(h^{2p})$ in the interior of the interface, and truncation error $\mathcal{O}(h^{p})$ on a few grid points near the edge of the interface. Therefore, in the semi–discretization (14), the largest truncation error is $\mathcal{O}(h^{p-2})$ introduced by the first three penalty terms in (15) and (16) because of $\tau,\boldsymbol{H_{ux}^{-1}},\boldsymbol{H_{uy}^{-1}}\sim\mathcal{O}(h^{-1})$. The truncation error $\mathcal{O}(h^{p-2})$ is only localized on a few number of grid points at the corner of two adjacent blocks. According to the accuracy analysis in [35], we may expect a rate of convergence $p+1$ of the semi–discretization (14), the same as the scheme developed in [36]. Note that the $p+1$ convergence rate is one order lower than the expected convergence rate when the grid interface is conforming. If the interpolation error at the edge could be improved to $\mathcal{O}(h^{p+1})$, then the expected rate of convergence would be $p+2$. However, it is proved in [20] that such norm–compatible interpolation operators cannot be constructed.

4.2 Numerical interface treatment of non–conforming blocks

An example of non–conforming blocks is shown in Figure 1b. The lower left corner of the upper right domain sits in the middle of the right boundary of the left domain. Such an interface configuration is sometimes called a T–junction interface.

For a T–junction interface, the interface conditions must be imposed on the glue grid, which is different from what is usually done for conforming blocks. The technique and the corresponding interpolation operators are constructed in [16], and are used for a T–junction interface on a Cartesian grid in [36]. On a curvilinear grid, the discretization is performed in a similar way, but the grid transformation must be done carefully.

Coordinate transformation is performed block–wise, therefore an interface between two blocks is transformed twice. A common strategy is to transform each block in the physical domain to the unit square in the reference domain. This works well with conforming blocks. However, with nonconforming blocks such as in Figure 1b the interfaces are transformed differently in different blocks. As a consequence, the transformed equation must be scaled as explained in [16] to obtain an energy stable scheme. The energy stable semi–discretization can then be constructed in a similar way as in the preceding sections.

5 Numerical experiments

In this section, numerical experiments are performed to verify the stability and accuracy property of the numerical schemes developed in this paper. The diagonal norm SBP operators used in the numerical experiments can be found in [31] for $D_{1}\approx\partial/\partial x$ and in [21] for $D_{2}^{(b)}\approx\partial/\partial x(b(x)\partial/\partial x)$. The interpolation operators can be found in [16, 22]. The L2 errors are computed as the norm of the difference between the exact solution $u_{ex}$ and the numerical solution $u_{h}$ according to

[TABLE]

where $h_{x}$ and $h_{y}$ are the mesh size in the $x$ and $y$  spatial direction, respectively.

5.1 An extreme interface

We consider the wave equation on the domain $[-1,1]\times[0,1]$, separated by the interface $x=4\sin(7\pi y)/5$. The domain and mesh are depicted in Figure 2a and 2b, respectively. The aim of this experiment is to verify that the scheme is stable even when an interface with a large curvature is present in the domain, but not to test accuracy or convergence rate. As can be seen in Figure 2b, the mesh is of bad quality due to large distortion.

The wave equation is discretized by the sixth order SBP operators in each subdomain, and patched together by the SAT method using the sixth order interpolation operators [22]. The semi–discretization can be written as a system of ordinary differential equations

[TABLE]

where $w$ is the numerical solution, $\boldsymbol{D}$ is the spatial discretization operator including boundary and interface terms, and $F$ corresponds to the forcing function and boundary data evaluated on the grid.

First, we use $21\times 21$ grid points in the left domain and $41\times 41$ grid points in the right domain, and perform an eigenvalue analysis. Stability requires that all the eigenvalues of $\boldsymbol{D}$ are real and non–positive. In Figure 3a, we plot the eigenvalues of $\boldsymbol{D}$ multiplied by the square of the mesh size in the right domain, denoted by $\lambda$, and observe that they are indeed real and non–positive.

Next, we test the scheme with a much finer mesh, $101\times 101$ grid points in the left domain and $201\times 201$ grid points in the right domain. Instead of an eigenvalue analysis, we perform a long time simulation by using the manufactured solution

[TABLE]

for initial and Neumann boundary data. We choose the classical Runge–Kutta method as the time integrator, and let the wave propagate for ten temporal periods. The L2 error at each time step is plotted in Figure 3b. We observe that the L2 error is bounded in time.

5.2 A T–junction interface

We consider the same domain $[-1,1]\times[0,1]$ as in the previous experiment, but with interfaces depicted in Figure 1b. The interface in the vertical direction is defined by $x=\sin(3\pi y/2)/3$. The intersection point $(\bar{x},\bar{y})$ of the two interfaces is chosen by letting $\bar{y}=0.621$. The interface in the horizontal direction is defined by $y=\sin(\pi x/2)/5+\bar{y}-\sin(\pi\sin(3\pi\bar{y}/2)/6)/5$. The numbers of grid points in the left, lower right and upper right domain are $26\times 52,26\times 26$ and $51\times 26$, respectively. Both the blocks and interfaces are non–conforming, see a close–up in Figure 4a. When refining the mesh, the number of grid points is doubled in each spatial direction in each domain.

To test accuracy and rate of convergence, we use the manufactured solution

[TABLE]

to obtain initial and Neumann boundary data, and propagate the wave until $t=2$. With this analytical solution, there is no forcing term in the equation.

We solve the equation by the fourth and sixth order SBP–SAT finite difference method, and use the classical Runge–Kutta method to integrate in time. Since the interface in this experiment is a T–junction, we use the fourth and sixth order accurate interpolation operators constructed in [16]. The time step is chosen small enough so that the error in the solution is determined by the spatial discretization. The errors in L2 norm are shown in Figure 4b, and the associated rates of convergence are given at the end of each error plot. In the figure, we use new SAT as the legend to denote the results obtained by the scheme in this paper, and old SAT to denote the result obtained by the scheme in [36]. The $x$–axis label $N$ is the number of grid points in the $x$–direction in the left domain.

We observe that the fourth and sixth order accurate scheme lead to third and fourth order convergence rate, respectively. This agrees well with the accuracy discussion in the end of Section 4.1 in this paper. We note that the sixth order method gives much smaller error than the fourth order methods. In addition, the four order method developed in this paper gives a smaller error than the fourth order method in [36]. We have also performed an experiment with the sixth order method in [36] and the numerical solution quickly blows up, indicating that method is unstable. This is not surprising because the energy analysis in [36] requires the norm–compatible condition of the interface operators, which are not satisfied by the sixth order operators.

6 Conclusion

We use the SBP–SAT finite difference method to solve the wave equation on a composite domain. The domain is divided by curved interfaces resulting in non–conforming blocks, and the grid is constructed in each block independently resulting in non–conforming grid interfaces. We develop new penalty terms to patch the blocks together by the SAT method. This extends the provably stable scheme from fourth order accuracy [36] to sixth order accuracy. Numerical experiments demonstrate the superiority of the new sixth order accurate scheme. In addition, we find that the new fourth order accurate scheme is more accurate than the fourth order accurate scheme in [36].

We note that eighth order and tenth order interpolation operators are constructed in [16], and can potentially be incorporated into the developed scheme in this paper. However, higher than sixth order accurate SBP operators for second derivative with variable coefficient have not yet been constructed.

Acknowledgement

The author would like to thank Professor Gunilla Kreiss at Uppsala University for the fruitful discussion and support.

Appendix

Stability is proved by the energy method, starting with multiplying the two semi–discretized equations in (14) from the left by $u_{t}^{T}\boldsymbol{H_{u}}$ and $v_{t}^{T}\boldsymbol{H_{v}}$, respectively, where $\boldsymbol{H_{u}}:=H_{ux}\otimes H_{uy}$ and $\boldsymbol{H_{v}}:=H_{vx}\otimes H_{vy}$. In the following derivation, we focus on the energy contribution from the first equation in (14) as the energy contribution from the second equation in (14) can be computed in a similar way.

First, we consider the energy contribution from the penalty term $SAT_{u}$:

[TABLE]

In the final energy estimate, we expect to see that the discrete energy is conserved in time. We note that the first part of the third term in (Appendix) can be written as the time derivative of a quadratic term

[TABLE]

With a positive $\tau$, the above term contributes positively to the discrete energy in terms of $\boldsymbol{E_{ux}}u$. By using the norm–compatible property of the interpolation operators, we find that the first part of the second term in (Appendix) can also be written as the time derivative of a quadratic term

[TABLE]

which contributes to the energy positively in terms of $(E_{ux}\otimes I_{u2v})u$. We observe that in (Appendix), there are also terms $\boldsymbol{E_{ux}S_{ux}}u$ and $\boldsymbol{E_{ux}D_{1uy}}u$. Therefore, we need the corresponding positive energy contributions, which come from the SBP operators as shown below.

The energy contribution from the two mixed–derivative terms is

[TABLE]

Note that we have used the equality (1) for $\boldsymbol{D_{1ux}}$ and $\boldsymbol{D_{1uy}}$, and have excluded boundary terms that do not correspond to the interface. On the right–hand side of (20), the first two terms are volume terms, involving the numerical solution $u$ in the entire domain $\Omega_{u}$; the third term is an interface term, involving $u$ on the interface.

The energy contribution from the SBP approximation is

[TABLE]

where the two remainder terms are $\overline{R^{(a)}_{ux}}=\sum_{i=1}^{n_{uy}}{R^{(a_{i})}_{ux}}\otimes E^{i}_{uy}$ and $\overline{R^{(c)}_{uy}}=\sum_{i=1}^{n_{ux}}{R^{(c_{i})}_{uy}}\otimes E^{i}_{ux}$. On the right–hand side of (Appendix), the third is an interface term, while the others are volume terms. For the volume terms in (20) and (Appendix), we have

[TABLE]

Next, we split $\begin{bmatrix}{\Lambda_{a}}&{\Lambda_{b}}\\ {\Lambda_{b}}&{\Lambda_{c}}\end{bmatrix}$ into two parts, and obtain

[TABLE]

where $\boldsymbol{I}$ is an identity operator, and $\delta$ is defined in (17). The first term in (22) is the time derivative of a non–positive quantity. From the second term in (22), we need to get quadratic terms for both $\boldsymbol{E_{ux}D_{1ux}}u$ and $\boldsymbol{E_{ux}S_{ux}}u$ on the interface. We write

[TABLE]

Note that in the above derivation, we use Lemma 3 to obtain (23), and Lemma 2 to obtain (24). We have obtained both the time derivative of quadratic terms for $\boldsymbol{E_{ux}S_{ux}}u$ and $\boldsymbol{E_{ux}D_{1ux}}u$.

After a very similar derivation of energy contribution for the second equation in (14), we move all terms to one side and write it in the form $\frac{d}{dt}\boldsymbol{G}=0$. The final step of the energy analysis is to determine the penalty parameter $\tau$ so that $\boldsymbol{G}$ is a discrete energy satisfying $\boldsymbol{G}\geq 0$. With some algebraic calculations, we may write the energy contribution from all interface terms as $x_{I}^{T}\boldsymbol{A}x_{I}$, where $x_{I}$ is the vector in the form

[TABLE]

The matrix $\boldsymbol{A}$ is a block matrix in the form $\boldsymbol{A}=\begin{bmatrix}\boldsymbol{A_{1}},&\boldsymbol{0},\\ \boldsymbol{0},&\boldsymbol{A_{2}}\end{bmatrix}$, where

[TABLE]

[TABLE]

We note that $\boldsymbol{A_{1}}$ is symmetric, and each submatrix of $\boldsymbol{A_{1}}$ is a diagonal matrix of dimension $n_{uy}$. Therefore, we can write $\boldsymbol{A_{1}}$ as a sum of $n_{uy}$ matrices, where the $i^{th}$ matrix takes the $i^{th}$ diagonal element of each submatrix with all the other elements zero. We then compute the eigenvalues of each matrix, and require them to be non–negative. By considering $\boldsymbol{A_{2}}$ in the same way, we obtain the limit on the penalty parameter for which $\boldsymbol{A}$ is positive semi–definite

[TABLE]

where $a_{\max}$, $b_{\max}$, $\alpha_{\max}$ and $\beta_{\max}$ are the maximum of functions $a(x,y)$, $b(x,y)$, $\alpha(x,y)$ and $\beta(x,y)$ evaluated on the interface, respectively.