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

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.


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-byparts 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 boxshaped 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 nonconforming. 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 normcontracting. 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 normcontracting. 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.

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, · · · , 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 ), · · · , U(x N )] T .

Definitions of SBP operators
The SBP concept and the first derivative SBP operator D 1 ≈ ∂/∂x are introduced in [19], and later refined in [31]. Formally it is defined as follows.
The operator H defines the SBP norm, and leads to the identity which is the discrete analogue of the integration-by-parts formula since the norm H is also a quadrature [7,14]. For the second derivative, we distinguish between a constant coefficient operator Definition 2. A difference operator D 2 = H −1 (−M + BS) approximating ∂ 2 /∂x 2 is a diagonal norm second derivative SBP operator if H is diagonal positive definite, M is symmetric positive semi-definite, B = diag(−1, 0, · · · , 0, 1), and the first and last row of S approximate ∂/∂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 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 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 , and the first and last row of S approximate ∂/∂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 where M is also symmetric positive semi-definite, σ > 0 is a constant independent of h, B and S are the same as in Definition 2, and 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 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 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 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 (b) 2 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 O(h 2p ) or O(h p ). In the discussion of accuracy, we make the truncation error of the scheme precise.

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

The wave equation on a Cartesian grid
We start by considering the wave equation in two space dimensions in a composite domain Ω = [0, 1] 2 with an interface Γ at x = 0.5. The left and right domain are denoted by Ω u and Ω v , respectively, and the equations are At the interface the physical conditions are 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 × n uy grid points in Ω u and n vx × n vy grid points in Ω v . We are particularly interested in a non-conforming interface when n uy = 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 Ω u to Ω v , and from Ω v to Ω u , respectively. In the SBP-SAT finite difference framework, these operators must satisfy certain conditions so that the scheme could be energy stable.
In [36], it is also defined that the interpolation operators are normcontracting if the two operators 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 O(h 2p ) in the interior of the interface, and 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 normcompatible 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 where and 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 Ω u and Ω v , respectively. Most operators in two space dimensions can be extended from the corresponding one dimensional operators by using a Kronecker product ⊗. 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 H −1 ux equals to H −1 ux ⊗ I uy , where H −1 ux is the inverse of the SBP norm in the x-direction acting on u, and I uy is an identity operator. The operator E extracts the numerical solution at the interface. D u u and 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 E ux u used in [33,36] is broken into two parts: 1/2E ux u in (8b) and 1/2E ux ⊗ I v2u I u2v u in (8c). Note the relation between them: E ux ⊗ I v2u I u2v u is just E ux u interpolated to the grid on the interface of Ω v , then interpolated back to the grid of Ω u . Since the interpolation is not exact, E ux ⊗I v2u I u2v u differs from 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 τ such that where θ is the constant in Lemma 1, and h ux and h vx are the mesh size in the x-direction in Ω u and Ω v , respectively.
Proof. We prove stability by the energy method. Multiplying from the left of (6) by u T t (H ux ⊗ H uy ) and (7) by v T t (H vx ⊗ H vy ), we obtain 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 τ 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 G as where Clearly, G 1 ≥ 0. By Young's inequality, we have G 2 ≥ 0 and G 3 ≥ 0 if τ ≥ 1/(2θh ux ) and τ ≥ 1/(2θ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 G in (11), with the inequality resulted from the norm-contracting condition.

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.

Numerical interface treatment of conforming blocks
With only conforming blocks in the domain, the corners of adjacent blocks meet. We consider again the domain Ω = [0, 1] 2 but partitioned into two blocks Ω u and Ω v by a curved interface. The grids are then constructed independently in each block, see an illustration in Figure 1a. The grids 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 where the Jacobians J U , J V > 0 and the variable coefficients satisfy ac−b 2 > 0 and αγ − β 2 > 0. The interface conditions become on (x, y) ∈ Γ. 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 a b b c and α β β γ because of ac − b 2 > 0 and αγ − β 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 where and We now clarify the notations in the semi-discretization (14).

The mixed-derivative terms: (bU
where Λ b is a diagonal matrix with diagonal entries b(x, y) evaluated on the grid. The operators approximating the other three mixedderivative terms are constructed in a similar way.
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 where a i is a(x i , y) evaluated on the grid, and E i uy has value one in entry (i, i) and zeros elsewhere. The second derivative operator D (a i ) 2ux 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.
Let a max , b max , α max and β max denote the maximum values of the variable coefficients a(x, y), b(x, y), α(x, y) and β(x, y) evaluated on the interface, respectively. We also denote h ux and h vx the mesh size in the x-direction in Ω u and Ω 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 τ satisfies where σ is defined in Lemma 2, and δ 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 τ 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 τ 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 O(h 2p ) in the interior and O(h p ) near the boundaries and the interface. The interpolation operators constructed in [16,22] have truncation error O(h 2p ) in the interior of the interface, and truncation error 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 O(h p−2 ) introduced by the first three penalty terms in (15) and (16) 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 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.

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.

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 ≈ ∂/∂x and in [21] for D The interpolation operators can be found in [16,22]. The L 2 errors are computed as the norm of the difference between the exact solution u ex and the numerical solution u h according to where h x and h y are the mesh size in the x and y spatial direction, respectively.  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.

An extreme interface
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 where w is the numerical solution, 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 × 21 grid points in the left domain and 41 × 41 grid points in the right domain, and perform an eigenvalue analysis. Stability requires that all the eigenvalues of D are real and non-positive. In Figure  3a, we plot the eigenvalues of D multiplied by the square of the mesh size in 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 L 2 error at each time step is plotted in Figure 3b. We observe that the L 2 error is bounded in time. To test accuracy and rate of convergence, we use the manufactured solution U = cos(3πx + 1) cos(4πy + 2) cos(5πt + 3).  Figure 3: (a) Eigenvalues of the spatial discretization operator (b) L 2 error in ten temporal periods 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 L 2 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.

A T-junction interface
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.

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. on the interface. We write