Non-conforming Interface Conditions for the Second-Order Wave Equation

Abstract

Imposition methods of interface conditions for the second-order wave equation with non-conforming grids is considered. The spatial discretization is based on high order finite differences with summation-by-parts properties. Previously presented solution methods for this problem, based on the simultaneous approximation term (SAT) method, have shown to introduce significant stiffness. This can lead to highly inefficient schemes. Here, two new methods of imposing the interface conditions to avoid the stiffness problems are presented: (1) a projection method and (2) a hybrid between the projection method and the SAT method. Numerical experiments are performed using traditional and order-preserving interpolation operators. Both of the novel methods retain the accuracy and convergence behavior of the previously developed SAT method but are significantly less stiff.

1 Introduction

It is well known that high order finite differences are highly efficient for large-scale wave propagation problems [1]. However, the design of such schemes requires particular care at the boundaries to obtain stability. One way to obtain stable high order finite difference schemes is to use finite difference operators with a summation-by-parts (SBP) property together with simultaneous-approximation-terms (SBP-SAT) [2], the projection method (SBP-P) [3, 4] or ghost points (SBP-GP) [5]. SBP finite difference operators are essentially standard finite difference stencils in the interior with boundary closures carefully designed to mimic integration by parts in the discrete setting. The SBP difference operators have an associated discrete inner product such that a discrete energy equation that is analogous to the continuous equation can be derived. The boundary conditions should be imposed such that the scheme exhibits no non-physical energy growth, sometimes referred to as strict stability [6]. The SAT method achieves this by adding penalty terms that weakly impose the boundary conditions such that the resulting scheme is stable, see for example [7]. The SBP-GP method adds ghost points at the boundaries and computes their values such that the boundary conditions are imposed and the scheme is stable [8, 9]. The projection method derives an orthogonal projection and rewrites the problem such that it is solved in the subspace of solutions where the boundary conditions are exactly fulfilled, see [10].

An important aspect of finite difference methods is the ability to split the computational domain into blocks and couple them across the interfaces. This is necessary to handle complex geometries, but also to increase the efficiency of the schemes. For example, in the case of the wave equation it may happen that a small wave speed in certain regions of the domain demands a high grid resolution. By decomposing the domain into blocks a finer grid spacing can be used only in these regions and a coarser grid where the wave speed is higher. This way the required number of grid points per wave length can be obtained, without using a unnecessarily refined grid in regions where its not needed. In general, the grid points at each side of an interface are non-conforming, in which case interpolations are used to couple the solutions. In the framework of SBP finite differences, it is crucial that the method of imposing the interpolated interface conditions preserves the SBP properties of the difference operators.

The construction of interpolation operators along with SATs to obtain stable schemes with non-conforming interfaces has received significant attention in the past [11,12,13]. In [11] so-called SBP-preserving interpolation operators (here referred to as norm-compatible) were first constructed and used to derive stable schemes for general hyperbolic and parabolic problems. However, it was noted in [13, 14] that the global convergence rate was decreased by one (compared to the convergence rate with conforming grids) for problems involving second derivatives in space. In [15] this is solved by constructing order-preserving (OP) interpolation operators along with SATs such that the global convergence rate is preserved. The new operators come in two norm-compatible pairs (a pair consists of two operators interpolating in opposite directions), where one of the operators in each pair is of one order higher accuracy at the boundaries. In Sect. 2.2 a more detailed discussion on the accuracy of the interpolation operators is presented.

A major downside of the SBP-SAT discretizations is the necessary decomposition of the second derivative SBP operator, often referred to as the “borrowing trick” [16], necessary to obtain an energy estimate. This procedure is known to introduce additional stiffness to the problem. The main contribution of the current work is two new methods avoiding this problem, one using SBP-P and the other a hybrid SBP-P-SAT. The analysis and numerical experiments are done on the second-order wave equation. However, the discrete Laplace operators presented are equally applicable to the heat equation and the Schrüdinger equation.

The paper is structured as follows: In Sect. 2 some necessary definitions and the discrete operators are introduced. In Sect. 3 the continuous problem is presented. The new semi-discrete schemes are presented in Sect. 4. The time discretization is presented in Sect. 5. In Sect. 6 numerical experiments validating the new methods and comparing them to the SBP-SAT schemes are presented. Conclusions are drawn in Sect. 7.

2 Definitions

Let

$$\begin{aligned} (f,g)_\Omega = \int _\Omega f g \, dx \quad \text {and} \quad \Vert f\Vert ^2_\Omega = (f,f), \end{aligned}$$

(1)

define the $L^2$-inner product and the corresponding norm for sufficiently smooth functions f, g on a two-dimensional domain $\Omega$. Let also

$$\begin{aligned} \langle f, g \rangle _{\partial \Omega } = \int _{\partial \Omega } f g \, dS \end{aligned}$$

(2)

denote integration along the boundary $\partial \Omega$. With this notation, the first Green’s identity reads

$$\begin{aligned} (f,\Delta g)_\Omega = \langle f, n \cdot \nabla g \rangle _{\partial \Omega } - (\nabla f, \nabla g)_\Omega , \end{aligned}$$

(3)

where n is the outward normal to the boundary.

2.1 Summation-by-Parts Finite Differences

We begin by considering one-dimensional SBP operators. Let $J = [x_l,x_r]$ define an interval and let ${\textbf{x}} = [x_1,x_2,...,x_m]$ be an equidistant discretization of J, where

$$\begin{aligned} x_i = x_l + (i-1) h, \quad i = 1,2,...,m, \quad \text {and} \quad h = \frac{x_r - x_l}{m-1}, \end{aligned}$$

(4)

is the step size. Let $f \in C^2(J)$ and let $f({\textbf{x}}) = [f(x_1),f(x_2),...,f(x_m)]$ be the restriction of f on ${\textbf{x}}$. Introduce the second derivative finite difference approximation matrix $D_2$ such that

$$\begin{aligned} D_2 f({\textbf{x}}) \approx \frac{d^2}{dx^2} f. \end{aligned}$$

(5)

The one-dimensional second derivative operator $D_2$ has summation-by-parts properties if it satisfies

$$\begin{aligned} D_2 = H^{-1}\left( -M + e_r d_r^\top - e_l d_l^\top \right) , \end{aligned}$$

(6)

where H is diagonal and positive definite, M is symmetric and positive semi-definite, $e_{l,r}^\top$ are row-vectors extracting the solution at the first and last grid points and $d_{l,r}^\top$ are row-vectors approximating the first derivative of the solution at the first and last grid points [17]. The matrix H defines a discrete inner product for vectors ${\textbf{f}},{\textbf{g}} \in {\mathbb {R}}^m$ as

$$\begin{aligned} \langle {\textbf{f}},{\textbf{g}}\rangle _H = {\textbf{f}}^\top H {\textbf{g}}. \end{aligned}$$

(7)

Fig. 1

Computational domain showing the orientation of the solution vectors

We now consider the discretization of the two-dimensional domain $\Omega$. Assume that $\Omega$ is square, and discretize it using m equidistant grid points in each dimension. In the present study general domains that would necessitate coordinate transforms are not considered (although this is highly relevant and an interesting extensions for a future study). Let ${\textbf{f}} = [{\textbf{f}}^{(1)},{\textbf{f}}^{(1)},...,{\textbf{f}}^{(m)}]^\top$, where ${\textbf{f}}^{(i)} = [f_1^{(i)},f_2^{(i)},...,f_m^{(i)}]$ contains the solution points along the y-axis at $x_i$, be a two-dimensional grid function in a column-major order, see Fig. 1. Let ${\textbf{g}}$ be a grid function constructed analogously. Note that now ${\textbf{f}},{\textbf{g}} \in {\mathbb {R}}^{m^2}$. The one-dimensional SBP operators are extended to two dimensions using Kronecker products as follows:

$$\begin{aligned} \begin{aligned} D_{2x}&= (D_2 \otimes I_{m}), \quad{} & {} D_{2y} = (I_{m} \otimes D_2), \\ H_x&= (H \otimes I_{m}),{} & {} H_y = (I_{m} \otimes H), \\ M_x&= (M \otimes I_{m}),{} & {} M_y = (I_{m} \otimes M), \\ e_W&= (e_l^\top \otimes I_{m}),{} & {} e_E = (e_r^\top \otimes I_{m}), \\ e_S&= (I_{m} \otimes e_l^\top ),{} & {} e_N = (I_{m} \otimes e_r^\top ), \\ d_W&= (d_l^\top \otimes I_{m}),{} & {} d_E = (d_r^\top \otimes I_{m}), \\ d_S&= (I_{m} \otimes d_l^\top ),{} & {} d_N = (I_{m} \otimes d_r^\top ), \\ \end{aligned} \end{aligned}$$

(8)

where $I_m$ is the $m \times m$ identity matrix. The subscripts in the first three rows of equations indicate the direction in which the operator is acting, for example $D_{2x}$ approximates the second derivative in the x-direction. The subscripts in the last four rows indicate the boundary, for example $e_W$ is used to extract the boundary values at the western boundary.

The discrete Laplace operator is given by

$$\begin{aligned} D_L = D_{2x} + D_{2y}. \end{aligned}$$

(9)

We also have the following two-dimensional discrete inner product and norm:

$$\begin{aligned} ({\textbf{f}},{\textbf{g}})_{{\bar{H}}} = {\textbf{f}}^\top {{\bar{H}}} {\textbf{g}} \quad \text {and} \quad \Vert {\textbf{f}}\Vert ^2_{{\bar{H}}} = ({\textbf{f}},{\textbf{f}})_{{\bar{H}}}, \end{aligned}$$

(10)

where ${\bar{H}} = H_x H_y$. Using the SBP properties of $D_2$ (6), the discrete two-dimensional Laplace operator can be written as

$$\begin{aligned} D_L = H_x^{-1} (-M_x + e_E^\top d_E - e_W^\top d_W) + H_y^{-1} (-M_y + e_N^\top d_N - e_S^\top d_S), \end{aligned}$$

(11)

or

$$\begin{aligned} \begin{aligned} ({\textbf{f}},D_L {\textbf{g}})_{{\bar{H}}} =&\, -{\textbf{f}}^\top (H_y M_x + H_x M_y) {\textbf{g}} + \langle e_E {\textbf{f}}, d_E {\textbf{g}}\rangle _H - \langle e_W {\textbf{f}}, d_W {\textbf{g}}\rangle _H \\&+\, \langle e_N {\textbf{f}}, d_N {\textbf{g}}\rangle _H - \langle e_S {\textbf{f}}, d_S {\textbf{g}}\rangle _H, \end{aligned} \end{aligned}$$

(12)

which is the discrete equivalent to Green’s identity (3) on $\Omega$.

The focus of the present work is on multiblock discretizations. Throughout the paper, a two-block domain $\Omega = \Omega _L \cup \Omega _R$, where $\Omega _L = [x_l,0]\times [y_l,y_r]$ is the left block and $\Omega _R = [0,x_r]\times [y_l,y_r]$ the right block. Here $x_{l,r}$ and $y_{l,r}$ define the outer boundary and an internal boundary, or interface, is located at $x = 0$. For notational simplicity, assume that both blocks are square, i.e., $x_r = - x_l = y_r - y_l$. Each block is discretized by an equidistant grid with $m^{(u)} \times m^{(u)}$ and $m^{(v)} \times m^{(v)}$ grid points in the left and right blocks, respectively. Throughout the paper, superscripts (u) and (v) will be used to denote operators in the left and right blocks. For example, $D_L^{(u)}$ is the discrete Laplace operator in the left block and $d_S^{(v)}$ approximates the y-derivative at the southern boundary in the right block. Solution vectors in the left and right blocks will be denoted ${\textbf{u}}$ and ${\textbf{v}}$, respectively.

For the SBP property (6) to hold, a decrease in order of accuracy of $D_2$ at a few grid points close to the boundaries is required. With the operators used here, for a 2pth-order accurate inner stencil, the order of accuracy at the boundaries is p. This will lead to a decrease in overall convergence rate. For one-dimensional problems it can be shown that the convergence rate is at least $\min (2p,p+2)$ for a point wise stable discretization of the second-order wave equation [18]. In [15], $\min (2p,p+2)$ is referred to as the ideal convergence rate for the two-dimensional problem with non-conforming interfaces. This is the rate that the OP interpolation operators are designed to preserve. In this paper, 4th and 6th order SBP operators and interpolation operators are used, i.e., $p = 2$ and $p = 3$. Thus, the convergence rates we hope to obtain are 4 and 5 with the OP interpolation operators, and 3 and 4 with traditional interpolation operators, for the 4th and 6th order SBP operators, respectively.

Remark 1

I have assumed that the two blocks, $\Omega _L$ and $\Omega _R$, are square and that each block is discretized using the same number of grid points in each dimension. This choice is made purely from a readability perspective. The theory presented is equally applicable to rectangular blocks discretized with different number of grid points in each dimension, but the analysis would become unnecessarily cluttered.

2.2 Interpolation Operators

If $m^{(u)} \ne m^{(v)}$, the grid points at the interface are non-conforming and interpolation operators are needed to couple the solutions between the two blocks. The numerical experiments presented in this paper is performed using a 1:2 grid ratio, i.e., the grid size in the left block is twice the grid size in the right block, see Fig. 2. Let $I_{u2v}$ denote the operator interpolating from left to right, and $I_{v2u}$ the operator interpolating from right to left. To obtain stability, we require that the pair of interpolation operators are norm-compatible, i.e., they must satisfy

$$\begin{aligned} (I_{v2u} {\textbf{v}}, {\textbf{u}})_{H^{(u)}} = ({\textbf{v}},I_{u2v} {\textbf{u}})_{H^{(v)}}, \quad \forall {\textbf{u}} \in {\mathbb {R}}^{m^{(u)}}, {\textbf{v}} \in {\mathbb {R}}^{m^{(v)}}. \end{aligned}$$

(13)

Interpolation operators satisfying (13) were first derived in [11], and first applied to the second-order wave equation in [13].

To satisfy the norm-compatibility condition (13), the order of accuracy of the interpolation operators close to the boundaries must be smaller than in the interior. As for the SBP operators, if the order of accuracy in the interior is 2p, at a few grid points close to the boundaries it will be p. In [13] it was indicated that this reduced order of accuracy close to the boundaries leads to a decrease of one in overall convergence rate when solving the wave equation on second-order form. An obvious remedy to this problem would be to increase the order of accuracy of the interpolation operators to $p+1$ close to the boundaries. But, as was shown in [19], this is not possible while also satisfying (13). In [15] this problem was avoided by noting that it is not necessary that both interpolation operators have order of accuracy $p+1$ close to the boundaries. To preserve the overall convergence rate, it is enough if one of them does and carefully chosen SATs are used. In [15], the new order-preserving interpolation operators and SAT schemes were used to prove stability for the wave equation on second order form, the heat equation, and Schrödingers equation, while simultaneously obtaining the full convergence rate $\min (2p,p+2)$. The OP interpolation operators come in two pairs: $I_{u2v}^{b}$ and $I_{v2u}^{g}$, and $I_{u2v}^{g}$ and $I_{v2u}^{b}$, where the ”good” operators (superscript g) have order of accuracy $p+1$ close to the boundaries and the ”bad” operators (superscript b) have order of accuracy p at a few points close to the boundaries. In the interior, all four operators have order of accuracy 2p. Each pair of the OP interpolation operators are norm-compatible, i.e., they satisfy (13).

Fig. 2

A two block domain with a 1:2 grid ratio non-conforming interface

3 Continuous Analysis

We consider the initial-value boundary problem

$$\begin{aligned} \begin{aligned} u_{tt}&= c_1^2 \Delta u, \quad{} & {} (x,y) \in \Omega _L,{} & {} {}&t \ge 0, \\ v_{tt}&= c_2^2 \Delta v,{} & {} (x,y) \in \Omega _R,{} & {} {}&t \ge 0, \\ n \cdot \nabla u&= g_u,{} & {} (x,y) \in \partial \Omega _L \setminus \partial \Omega _I, \quad{} & {} {}&t \ge 0, \\ n \cdot \nabla v&= g_v,{} & {} (x,y) \in \partial \Omega _R \setminus \partial \Omega _I, \quad{} & {} {}&t \ge 0, \\ u&= v,{} & {} (x,y) \in \partial \Omega _I,{} & {} {}&t \ge 0, \\ c_1^2 u_x&= c_2^2 v_x,{} & {} (x,y) \in \partial \Omega _I,{} & {} {}&t \ge 0, \\ \end{aligned} \end{aligned}$$

(14)

with initial data for u, $u_t$, v, and $v_t$ at $t = 0$. Here $\partial \Omega _{L,R}$ denote the boundaries of the blocks, $\partial \Omega _I$ denotes the interface, n is the outward pointing normal, $g_{u,v}$ are boundary data, and $c_1$ and $c_2$ are real constants.

Taking the $L_2$-inner product between $u_t$ and the first equation in (14), adding the result to the $L_2$-inner product between $v_t$ and the second equation in (14), and using (3) lead to the energy equation

$$\begin{aligned} \frac{d}{dt} E = 2c_1^2 \langle u_t, n \cdot \nabla u \rangle _{\Omega _L} + 2c_2^2 \langle v_t, n \cdot \nabla v \rangle _{\Omega _R}. \end{aligned}$$

(15)

The energy is given by

$$\begin{aligned} E = \Vert u_t\Vert ^2_{\Omega _L} + \Vert v_t\Vert ^2_{\Omega _R} + c_1^2 \Vert \nabla u\Vert ^2_{\Omega _L} + c_2^2 \Vert \nabla v\Vert ^2_{\Omega _R}, \end{aligned}$$

(16)

Inserting the interface and boundary conditions (the last four equations in (14)) and assuming $g_{u,v} = 0$ leads to energy conservation,

$$\begin{aligned} \frac{d}{dt} E = 0. \end{aligned}$$

(17)

This energy estimate is sufficient to show that (14) is stable and has a unique solution.

4 Spatial Discretization

We now turn to the spatial discretization, time is left continuous for now. Before discussing the discretization of (14), a short introduction to the projection method is included in Sect. 4.1. Then, in Sect. 4.2, the new discretizations of the multi-block problem (14) are presented.

4.1 The Projection Method

Consider the general linear constrained second-order ODE

$$\begin{aligned} \begin{aligned} {\textbf{w}}_{tt}&= D {\textbf{w}}, \quad{} & {} t> 0, \\ L {\textbf{w}}&= 0,{} & {} t > 0, \\ {\textbf{w}}&= \mathbf {f_1},{} & {} t = 0, \\ {\textbf{w}}_t&= \mathbf {f_2},{} & {} t = 0, \end{aligned} \end{aligned}$$

(18)

where ${\textbf{w}} \in {\mathbb {R}}^N$ and N is a positive integer. The matrix $D \in {\mathbb {R}}^{N \times N}$ may be the spatial discretization of some PDE and $L \in {\mathbb {R}}^{N \times s}$ the discretization of the boundary conditions (s equations). We wish to impose the constraint $L {\textbf{w}} = 0$ using the projection method. Let $V:= \left\{ {\textbf{w}} \in {\mathbb {R}}^N: L {\textbf{w}} = 0 \right\}$ be the subspace of vectors satisfying the constraint and define P as the projection from ${\mathbb {R}}^N$ onto V. The projection method imposes the constraint by solving the augmented ODE

$$\begin{aligned} \begin{aligned} {\textbf{w}}_{tt}&= P D P {\textbf{w}}, \quad{} & {} t > 0, \\ {\textbf{w}}&= \mathbf {f_1},{} & {} t = 0, \\ {\textbf{w}}_t&= \mathbf {f_2},{} & {} t = 0. \end{aligned} \end{aligned}$$

(19)

From (19) we immediately see that ${\textbf{w}}_{tt} = P {\textbf{w}}_{tt}$, and by integrating in time twice we obtain ${\textbf{w}} = P {\textbf{w}} \Leftrightarrow L {\textbf{w}} = 0$ if $\mathbf {f_{1,2}} \in V$. We conclude that if the initial data satisfy the constraint, so will the solution to (19).

The particular form of P is determined by considering the stability properties of (19). Assume that we have some inner product $(\cdot ,\cdot )_H$ associated with D and apply the energy method, we obtain

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert {\textbf{w}}_t\Vert ^2_H = ({\textbf{w}}_t,PDP {\textbf{w}})_H. \end{aligned}$$

(20)

If P is self-adjoint with respect to the inner product, i.e.,

$$\begin{aligned} ({\textbf{w}}_1,P{\textbf{w}}_2)_H = (P{\textbf{w}}_1,{\textbf{w}}_2)_H, \quad \forall {\textbf{w}}_{1,2} \in {\mathbb {R}}^N, \end{aligned}$$

(21)

we can rewrite (20) as

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert {\textbf{w}}_t\Vert ^2_H = (\hat{{\textbf{w}}}_t,D \hat{{\textbf{w}}})_H, \end{aligned}$$

(22)

where $\hat{{\textbf{w}}} = P {\textbf{w}}$. Obviously the constraint is fulfilled for $\hat{{\textbf{w}}}$, $L \hat{{\textbf{w}}} = 0$. And, if $\mathbf {f_{1,2}} \in V$, we have ${\textbf{w}} = \hat{{\textbf{w}}}$. By using the SBP properties of D and choosing L so that the terms corresponding to energy growth are cancelled, a non-growing energy estimate is obtained immediately from (22).

The self-adjoint condition (21) is equivalent to P being an orthogonal projection with respect to H. Based on this, the following explicit formula for P can be derived:

$$\begin{aligned} P = I - H^{-1} L^\top ( L H^{-1} L^\top )^{-1} L, \end{aligned}$$

(23)

where I is the $N \times N$ identity matrix.

The projection method has been applied to various problems in the past. For more details see [3, 4, 10, 20]. See also [21] for examples of the projection method used for interface conditions.

4.2 Multi-block Analysis

We now consider the multi-block problem (14). The Neumann boundary conditions are imposed using the SAT technique, which has been extensively used in the past, see for example [22]. To make the analysis more readable, the terms corresponding to outer boundaries are left out. In Appendix 1, the semi-discrete schemes including boundary treatments are specified.

Let

$$\begin{aligned} {\textbf{w}} = \begin{bmatrix} {\textbf{u}} \\ {\textbf{v}} \end{bmatrix}, \end{aligned}$$

(24)

be the global semi-discrete solution vector of size $N = m^{(u)}m^{(u)} + m^{(v)}m^{(v)}$. Discretizing (14) in space without considering the interface conditions yields

$$\begin{aligned} \begin{aligned} {\textbf{w}}_{tt}&= {\hat{D}}_L {\textbf{w}}, \quad{} & {} t > 0, \\ {\textbf{w}}&= \mathbf {f_1},{} & {} t = 0, \\ {\textbf{w}}_t&= \mathbf {f_2},{} & {} t = 0, \end{aligned} \end{aligned}$$

(25)

where

$$\begin{aligned} {\hat{D}}_L = \begin{bmatrix} c_1^2 D_L^{(u)} &{} 0 \\ 0 &{} c_2^2 D_L^{(v)}\end{bmatrix}, \end{aligned}$$

(26)

and $\mathbf {f_{1,2}}$ is the discrete initial data. The global inner product associated with D is given by

$$\begin{aligned} ({\textbf{w}}_1,{\textbf{w}}_2)_{{\hat{H}}} = {\textbf{w}}_1^\top {\hat{H}} {\textbf{w}}_2, \quad \forall {\textbf{w}}_{1,2} \in {\textbf{R}}^N, \end{aligned}$$

(27)

where

$$\begin{aligned} {\hat{H}} = \begin{bmatrix} \bar{H}^{(u)} &{} 0 \\ 0 &{} \bar{H}^{(v)} \end{bmatrix}, \end{aligned}$$

(28)

We now wish to discretize the continuous interface conditions,

$$\begin{aligned} u = v, \quad c_1^2 u_x = c_2^2 v_x, \quad (x,y) \in \partial \Omega _I, \quad t \ge 0, \end{aligned}$$

(29)

using interpolation operators, and modify (25) using either SBP-P or SBP-P-SAT so that the discrete conditions are imposed and the modified ODE system is stable.

4.3 Stability with SBP-P

We begin by considering only the projection method to impose the interface conditions. The following lemma is the first main result of this paper:

Lemma 1

The ODE (19) with $D = {\hat{D}}_L$, $H = {\hat{H}}$, and L given by

$$\begin{aligned} L = \left[ {\begin{array}{*{20}l}{e^{(u)}_E} &{} {-I_{v2u} e^{(v)}_W} \\ {c_1^2 I_{u2v} d_E^{(u)}} &{} {-c_2^2 d_W^{(v)}}\\ \end{array} } \right] , \end{aligned}$$

(30)

is a stable approximation of (14) with a non-conforming interface.

Proof

Let $\hat{{\textbf{w}}} = \begin{bmatrix} \hat{{\textbf{u}}} \\ \hat{{\textbf{v}}} \end{bmatrix} = P {\textbf{w}}$ denote the projected solution vector. We begin by restating (22) for this problem. We have

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert {\textbf{w}}_t\Vert ^2_{{\hat{H}}} = c_1^2 (\hat{{\textbf{u}}}_t,D_L^{(u)} \hat{{\textbf{u}}})_{{\bar{H}}^{(u)}} + c_2^2(\hat{{\textbf{v}}}_t,D_L^{(v)} \hat{{\textbf{v}}})_{{\bar{H}}^{(v)}}. \end{aligned}$$

(31)

Using the SBP properties of $D_L^{(u,v)}$ (12), and rearranging results in

$$\begin{aligned} \frac{d}{dt} E = 2 c_1^2 (e_E^{(u)} \hat{{\textbf{u}}}_t, d_E^{(u)} \hat{{\textbf{u}}})_{H^{(u)}} - 2 c_2^2 (e_W^{(v)} \hat{{\textbf{v}}}_t, d_W^{(v)} \hat{{\textbf{v}}})_{H^{(v)}}, \end{aligned}$$

(32)

where E is an energy given by

$$\begin{aligned} \begin{aligned} E =&\, \Vert {\textbf{u}}_t\Vert ^2_{{\bar{H}}^{(u)}} + c_1^2 \hat{{\textbf{u}}}^ \top (H_y^{(u)} M_x^{(u)} + H_x^{(u)} M_y^{(u)}) \hat{{\textbf{u}}} \\&+\, \Vert {\textbf{v}}_t\Vert ^2_{{\bar{H}}^{(v)}} + c_2^2 \hat{{\textbf{v}}}^\top (H_y^{(v)} M_x^{(v)} + H_x^{(v)} M_y^{(v)}) \hat{{\textbf{v}}} \ge 0. \end{aligned} \end{aligned}$$

(33)

Since $L \hat{{\textbf{w}}} = L P {\textbf{w}} = 0$ by construction, the discrete interface conditions

$$\begin{aligned} e^{(u)}_E \hat{{\textbf{u}}} = I_{v2u} e^{(v)}_W \hat{{\textbf{v}}} \quad \text {and} \quad c_2^2 d_W^{(v)} \hat{{\textbf{v}}} = c_1^2 I_{u2v} d_E^{(u)} \hat{{\textbf{u}}}, \end{aligned}$$

(34)

hold exactly. Substituted into (32) results in

$$\begin{aligned} \frac{d}{dt} E = 2 c_1^2 ((I_{v2u} e^{(v)}_W \hat{{\textbf{v}}}_t, d_E^{(u)} \hat{{\textbf{u}}})_{H^{(u)}} - (e^{(v)}_W \hat{{\textbf{v}}}_t, I_{u2v} d_E^{(u)} \hat{{\textbf{u}}})_{H^{(v)}}). \end{aligned}$$

(35)

Using that the interpolation operators are norm-compatible, i.e., that they satisfy (13), we get

$$\begin{aligned} \frac{d}{dt} E = 0, \end{aligned}$$

(36)

which proves stability. $\square$

Remark 2

The scheme (19) with $D = {\hat{D}}_L$ and L given by (30) is the same as in [21] if the grids at the interface are conforming, and the interpolation operators are replaced with identity matrices.

4.4 Stability with SBP-P-SAT

With the hybrid method, instead of projecting the second interface condition

$$\begin{aligned} c_1^2 I_{u2v} d_E^{(u)} {{\textbf{u}}} = c_2^2 d_W^{(v)} {{\textbf{v}}}, \end{aligned}$$

(37)

it is imposed weakly using the SAT method. Let

$$\begin{aligned} {\tilde{D}} = {\hat{D}}_L + \begin{bmatrix} 0 &{} 0 \\ -(H_x^{(v)})^{-1} e^{(v) \top }_W c_1^2 I_{u2v} d_E^{(u)} &{} (H_x^{(v)})^{-1} e^{(v) \top }_W c_2^2 d_W^{(v)} \end{bmatrix}, \end{aligned}$$

(38)

denote a modified spatial operator. The second term in (38) is a SAT imposing (37) weakly on the equation for ${\textbf{v}}$. The following lemma is the second main result of this paper:

Lemma 2

The ODE (19) with $D = {\tilde{D}}$, $H = {\hat{H}}$, and L given by

$$\begin{aligned} L = \begin{bmatrix} e^{(u)}_E&-I_{v2u} e^{(v)}_W \end{bmatrix}, \end{aligned}$$

(39)

is a stable approximation of (14) with a non-conforming interface.

Proof

Let $\hat{{\textbf{w}}} = \begin{bmatrix} \hat{{\textbf{u}}} \\ \hat{{\textbf{v}}} \end{bmatrix} = P {\textbf{w}}$ denote the projected solution vector. We begin by restating (22) for this problem. We have

$$\begin{aligned} \begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert {\textbf{w}}_t\Vert ^2_{{\hat{H}}} =&\, c_1^2 (\hat{{\textbf{u}}}_t,D_L^{(u)} \hat{{\textbf{u}}})_{{\bar{H}}^{(u)}} + c_2^2(\hat{{\textbf{v}}}_t,D_L^{(v)} \hat{{\textbf{v}}})_{{\bar{H}}^{(v)}} \\&+\, \langle e_W^{(v)} \hat{{\textbf{v}}}, c_2^2 d_W^{(v)} \hat{{\textbf{v}}} - c_1^2 I_{u2v} d_E^{(u)} \hat{{\textbf{u}}} \rangle _{H^{(v)}} \end{aligned} \end{aligned}$$

(40)

Using the SBP properties of $D_L^{(u,v)}$ (12), and rearranging results in

$$\begin{aligned} \frac{d}{dt} E = 2 c_1^2 \left( \left( e^{(u)}_E \hat{{\textbf{u}}}_t, d^{(u)}_E \hat{{\textbf{u}}}\right) _{H^{(u)}} - \left( e^{(v)}_W \hat{{\textbf{v}}}_t,I_{u2v} d^{(u)}_E \hat{{\textbf{u}}}\right) _{H^{(v)}}\right) , \end{aligned}$$

(41)

where E is an energy given by

$$\begin{aligned} \begin{aligned} E =&\, \Vert {\textbf{u}}_t\Vert ^2_{{\bar{H}}^{(u)}} + c_1^2 \hat{{\textbf{u}}}^\top (H_y^{(u)} M_x^{(u)} + H_x^{(u)} M_y^{(u)}) \hat{{\textbf{u}}} \\&+\, \Vert {\textbf{v}}_t\Vert ^2_{{\bar{H}}^{(v)}} + c_2^2 \hat{{\textbf{v}}}^\top (H_y^{(v)} M_x^{(v)} + H_x^{(v)} M_y^{(v)}) \hat{{\textbf{v}}} \ge 0. \end{aligned} \end{aligned}$$

(42)

Since $L \hat{{\textbf{w}}} = L P {\textbf{w}} = 0$ by construction, the discrete interface condition

$$\begin{aligned} e^{(u)}_E \hat{{\textbf{u}}} = I_{v2u} e^{(v)}_W \hat{{\textbf{v}}}, \end{aligned}$$

(43)

holds exactly. Substituted into (41) results in

$$\begin{aligned} \frac{d}{dt} E = 2 c_1^2 \left( \left( I_{v2u} e^{(v)}_W \hat{{\textbf{v}}}_t, d_E^{(u)} \hat{{\textbf{u}}}\right) _{H^{(u)}} - \left( e^{(v)}_W \hat{{\textbf{v}}}_t, I_{u2v} d_E^{(u)} \hat{{\textbf{u}}}\right) _{H^{(v)}}\right) . \end{aligned}$$

(44)

Using that the interpolation operators are norm-compatible, i.e., that they satisfy (13), we get

$$\begin{aligned} \frac{d}{dt} E = 0, \end{aligned}$$

(45)

which proves stability. $\square$

Remark 3

With both SBP-P and SBP-P-SAT the key to obtaining energy stability is to interpolate the two interface conditions in opposite directions. In Sects. 4.3 and 4.4, continuity of the solution is imposed by interpolating right to left and the continuity of the flux by interpolating left to right. Conservative energy estimates can also be obtained by swapping the interpolations and using the transpose of (13). With only SBP-P this amounts to

$$\begin{aligned} L = \begin{bmatrix} I_{u2v} e^{(u)}_E &{} -e^{(v)}_W \\ c_1^2 d_E^{(u)} &{} -c_2^2 I_{v2u} d_W^{(v)} \end{bmatrix}, \end{aligned}$$

(46)

and with SBP-P-SAT

$$\begin{aligned} L = \begin{bmatrix} I_{u2v} e^{(u)}_E&-e^{(v)}_W \end{bmatrix}, \end{aligned}$$

(47)

and

$$\begin{aligned} {\tilde{D}} = {\hat{D}}_L + \begin{bmatrix} -(H_x^{(u)})^{-1} e^{(u) \top }_E c_1^2 d^{(u)}_E &{} (H_x^{(u)})^{-1} e^{(u) \top }_E c_2^2 I_{v2u} d^{(v)}_W\\ 0 &{} 0 \end{bmatrix}. \end{aligned}$$

(48)

Numerical experiments show that the differences between the choices in terms of accuracy and stiffness are minor, and dependent on the SBP and interpolation operators used. The results presented in this paper are obtained using the discretizations presented in Sects. 4.3 and 4.4.

4.5 Order-Preserving Interpolation

In [15] the order-preserving scheme is obtained by interpolating continuity of the solution using the ”good” operators and continuity of the flux using the ”bad” operators. With this choice, it turns out that the local truncation errors of the particular SATs used are balanced, meaning that they have the same minimal scaling with respect to the grid size. A detailed analysis of the truncation errors with SBP-P and SBP-P-SAT is out of scope in the present work. But, in the numerical experiments I investigate whether the OP interpolation operators can be used with the new methods to obtain a higher overall convergence rate as well. For the discretizations in Sects. 4.3 and 4.4, this amounts to replacing $I_{v2u}$ with $I_{v2u}^g$ and $I_{u2v}$ with $I_{u2v}^b$ in (30), (38), and (39). Note that with SBP-P and SBP-P-SAT only one pair of the OP interpolation operators is used, whereas the SBP-SAT discretization requires both pairs.

5 Time Discretization

All methods considered can be written as a system of second-order ODEs, given by

$$\begin{aligned} \begin{aligned} {\textbf{w}}_{tt}&= Q {\textbf{w}} + G(t), \quad{} & {} t > 0,\\ {\textbf{w}}(t)&= {\textbf{f}}_1,{} & {} t = 0, \\ {\textbf{w}}_t(t)&= {\textbf{f}}_2,{} & {} t = 0, \\ \end{aligned} \end{aligned}$$

(49)

where Q is a matrix approximating the spatial derivatives including boundary and interface conditions and G(t) contains the boundary data. In Appendix 1Q and G(t) are specified for all schemes used to produce the numerical results. In this paper (49) is solved using an explicit 4th order time-marching scheme [21, 23], given by

$$\begin{aligned} \begin{aligned} {\textbf{w}}^{(0)} =&\, {\textbf{f}}_1, \\ {\textbf{w}}^{(1)} =&\, \left( I + \frac{k^2}{2} Q \right) {\textbf{f}}_1 + k \left( I + \frac{k^2}{6} Q \right) {\textbf{f}}_2 + \frac{k^2}{2} G(0) + \frac{k^3}{6}G_t(0), \\ {\textbf{w}}^{(n+1)} =&\, \left( 2I + k^2 Q + \frac{k^4}{12} Q^2 \right) {\textbf{w}}^{(n)} - {\textbf{w}}^{(n-1)} \\&+ k^2 \left( I + \frac{k^2}{12} Q \right) G(t_n) + \frac{k^4}{12} G_{tt}(t_n), \end{aligned} \end{aligned}$$

(50)

where I is the identity matrix, k denotes the time step, and $t_n = nk$, $n = 0,1,...$, is the discrete time-level. In [23] it is shown that the scheme is stable if

$$\begin{aligned} k^2 \rho (Q) < 12, \end{aligned}$$

(51)

where $\rho (Q)$ denotes the spectral radius of Q. Introducing the undivided matrix ${\tilde{Q}} = h^{(v) ^2} Q$, where $h^{(v)}$ is the spatial interval in the right block, we get the stability condition

$$\begin{aligned} k < \sqrt{\frac{12}{\rho ({\tilde{Q}})}} h^{(v)}. \end{aligned}$$

(52)

The scaled spectral radius $\rho ({\tilde{Q}})$ depends on the discretization method, but not on the spatial interval $h^{(v)}$ (for large enough problems). Therefore, comparing the scaled spectral radius of the methods gives a good indication of the required time steps, and consequently the overall efficiency of the schemes.

6 Numerical Experiments

In this section numerical experiments are presented comparing the new discretizations to the SAT discretizations presented in [14] (traditional interpolation) and [15] (OP interpolation). The Neumann boundary conditions are imposed using the SAT method. All schemes are presented in Appendix 1. The domain is given by $[-10,10] \times [0,10]$ with an interface at $x = 0$. The left and right blocks are discretized with m and $2m - 1$ grid points in each dimension, which results in a non-conforming interface with grid ratio 1:2 at $x = 0$. The methods are compared in terms of efficiency (measured by the spectral radius) in Sect. 6.1, and accuracy for a problem with a known analytical solution in Sect. 6.2.

6.1 Spectral Radius

To measure the influence of material parameter jumps on the required time step the scaled spectral radius $\rho ({\hat{Q}})$ is plotted against C in Fig. 3, where $c_1 = 1$ and $c_2 = C$. Also plotted is the scaled spectral radius for a single block problem without an interface with wave speed $c = \max (C,0.5)$. Note that the wave speed for the single block problem is bounded from below by 0.5 since for $C < 0.5$ the scaled spectral radius of the two-block problem is dominated by the left block discretization. And hence we cannot hope to obtain a lower spectral radius than with the single block problem and $C = 0.5$. The results in Fig. 3 show that for all operators considered (4th and 6th order accurate with traditional and OP interpolation operators) the scaled spectral radius with SBP-SAT is significantly higher than with SBP-P and SBP-P-SAT. Note also that the scaled spectral radius with SBP-P and SBP-P-SAT exactly matches the single block discretization, indicating that the spectral radii with SBP-P and SBP-P-SAT are unaffected by the interface coupling procedure. In Table 1 the scaled spectral radii of the schemes are presented for $C = 0.5$, which with a 1:2 grid ratio corresponds to an equal number of grid points per wave length in each block. With this choice, for the 6th order OP operators with a given grid resolution, an approximately 2.5 times larger time step can be used with SBP-P or SBP-P-SAT compared to SBP-SAT. For the 4th order OP operators, the ratio is approximately 6.6.

Remark 4

The SBP-SAT schemes involve tuning the value of a parameter. Typically, increasing its value leads to a more accurate scheme (up to a point) at the cost of increasing the spectral radius. How to choose this parameter is not obvious, and one unclear aspect of the SAT method. The results in this paper are obtained using the same values as in [14, 15], see Appendix 1.

Table 1 Scaled spectral radius of RHS matrix with order-preserving and traditional interpolation operators with projection (SBP-P), hybrid projection and SAT (SBP-P-SAT), and SAT (SBP-SAT) for 4th and 6th order SBP operators

Fig. 3

Scaled spectral radius $\rho ({\hat{Q}})$ plotted versus C (ratio between wave speeds in fine block and coarse block) with order-preserving and traditional interpolation operators discretized using projection (SBP-P), hybrid projection and SAT (SBP-P-SAT), and SAT (SBP-SAT) for the 4th and 6th order SBP operators

6.2 Accuracy

In this section the accuracy of the methods is compared using an analytical solution given by

$$\begin{aligned} \begin{aligned} u&= \cos (x + y - \sqrt{2} c_1 t) + k_2 \cos (x - y + \sqrt{2} c_1 t), \\ v&= (1 + k_2) \cos (k_1 x + y {-} \sqrt{2} c_1 t), \end{aligned} \end{aligned}$$

(53)

where $k_1 = \sqrt{2 c_1^2/c_2^2-1}$ and $k_2 = (c_1^2 - c_2^2 k_1)/(c_1^2 + c_2^2 k_1)$. The wave speeds are set to $c_1 = 1$ and $c_2 = 0.5$. The boundary and initial data are given by (53). The time step is chosen as one tenth of the largest stable time step (with this choice the temporal errors are insignificant in comparison to the spatial errors). The convergence rate is approximated as

$$\begin{aligned} q = \frac{\log \left( \frac{e_1}{e_2}\right) }{\log \left( \frac{m_1}{m_2}\right) }, \end{aligned}$$

(54)

where $e_1$ and $e_2$ are errors in the H-norm (28) at $t = 2$ of two separate simulations with $m = m_1$ and $m = m_2$.

In Table 2 the error and convergence results of the SBP-P, SBP-P-SAT, and SBP-SAT discretizations are presented for the 4th and 6th order traditional and OP interpolation operators. Overall the accuracy of the SBP-P, SBP-P-SAT, and SBP-SAT schemes are very similar. With the traditional interpolation operators, 3rd and 4th order convergence are obtained with the 4th and 6th order operators respectively. And, with the order-preserving interpolation operators, convergence rates 4 and 5 are obtained. This shows that SBP-P and SBP-P-SAT exhibit the same convergence behaviors as previously observed with SBP-SAT, where the traditional interpolation operators lead to an order reduction whereas the OP interpolation operators retain the full convergence rates. One stand-out result is the accuracy with the 6th order traditional interpolation operators. With SBP-P and SBP-P-SAT, the errors with $m = 801$ are smaller by almost one magnitude compared to the error with SBP-SAT.

Table 2 Error (in base 10 logarithm) and convergence of 4th and 6th order traditional and order-preserving interpolation and SBP operators with projection (subscript p), hybrid (subscript h), and SAT (subscript s) discretizations

7 Conclusion

Two new SBP finite difference discretizations of the second-order wave equation with non-conforming grid interfaces are presented. The first scheme utilizes the projection method to impose the interface conditions and the second scheme a hybrid projection-SAT method. Semi-discrete energy conservation is shown for both discretizations using the energy method. Numerical experiments with traditional and order-preserving interpolation operators demonstrate similar accuracy and convergence behavior as for the SAT schemes. The most significant advantage of the new methods compared to SAT is the reduced spectral radius of the spatial operators. The new methods are less stiff than the SAT schemes, allowing for several times larger time steps with explicit time integration methods. Furthermore, it is found that the stiffness of the new schemes is the same as without the interface altogether, i.e., it is unaffected by the coupling procedure. Although the analysis and numerical experiments are done for the second-order wave equation, the discrete Laplace operator presented here can be directly applied to the heat equation and the Schrödinger equation. In a future study, the ideas introduced in this paper will be extended to general hyperbolic systems.

A Discretizations

The spatial discretization matrix Q and data vector G(t) in (49) are here presented for the single-block discretization with SBP-SAT and the two-block discretizations using SBP-P, SBP-P-SAT, and SBP-SAT.

1.1 A1. Single Block SBP-SAT

The single-block discretization with Neumann boundary conditions imposed using SBP-SAT is given by

$$\begin{aligned} Q = D_L + {\bar{\tau }}_W d_W + {\bar{\tau }}_E d_E + {\bar{\tau }}_S d_S + \bar{\tau }_N d_N, \end{aligned}$$

(55)

and

$$\begin{aligned} G(t) = -{\bar{\tau }}_W {\textbf{g}}_W(t) - {\bar{\tau }}_E {\textbf{g}}_E - {\bar{\tau }}_S {\textbf{g}}_S - {\bar{\tau }}_N {\textbf{g}}_N. \end{aligned}$$

(56)

Here

$$\begin{aligned} \begin{aligned} {\bar{\tau }}_W&= c^2 {\bar{H}}_x^{-1} e_W^\top , \quad {\bar{\tau }}_E = -c^2 {\bar{H}}_x^{-1} e_E^\top \\ {\bar{\tau }}_S&= c^2 {\bar{H}}_y^{-1} e_S^\top , \quad {\bar{\tau }}_N = -c^2 {\bar{H}}_y^{-1} e_N^\top , \end{aligned} \end{aligned}$$

(57)

and ${\textbf{g}}_{W,E,S,N}$ are vectors of g evaluated on the boundary grid points.

1.2 A2. Two-Block SBP-P

The two-block discretization with Neumann boundary conditions imposed using SBP-SAT and interface conditions imposed using SBP-P is given by

$$\begin{aligned} Q = P \left( \begin{bmatrix} c_1^2 D_L^{(u)} &{} 0 \\ 0 &{} c_2^2 D_L^{(v)} \end{bmatrix} + SAT_{BC} \right) P, \end{aligned}$$

(58)

and

$$\begin{aligned} G(t) = -P \begin{bmatrix} {\bar{\tau }}_W^{(u)} {\textbf{g}}_W^{(u)}(t) + {\bar{\tau }}_S^{(u)} {\textbf{g}}_S^{(u)}(t) + {\bar{\tau }}_N^{(u)} {\textbf{g}}_N^{(u)}(t) &{} 0 \\ 0 &{} {\bar{\tau }}_E^{(v)} {\textbf{g}}_E^{(v)}(t) + {\bar{\tau }}_S^{(v)} {\textbf{g}}_S^{(v)}(t) + {\bar{\tau }}_N^{(v)} {\textbf{g}}_N^{(v)}(t) \end{bmatrix}. \end{aligned}$$

(59)

Here

$$\begin{aligned} SAT_{BC} = \begin{bmatrix} {\bar{\tau }}_W^{(u)} d_W^{(u)} + \bar{\tau }_S^{(u)} d_S^{(u)} + {\bar{\tau }}_N^{(u)} d_N^{(u)} &{} 0 \\ 0 &{} \bar{\tau }_E^{(v)} d_E^{(v)} + {\bar{\tau }}_S^{(v)} d_S^{(v)} + \bar{\tau }_N^{(v)} d_N^{(v)} \end{bmatrix}, \end{aligned}$$

(60)

where the penalty matrices ${\bar{\tau }}^{(u)}_{W,S,N}$ and ${\bar{\tau }}^{(v)}_{E,S,N}$ are given by (57) evaluated in each respective block and ${\textbf{g}}^{(u)}_{W,S,N}$ and ${\textbf{g}}^{(v)}_{E,S,N}$ are boundary data. The projection operator is given by

$$\begin{aligned} P = I - {\hat{H}}^{-1} L^\top (L {\hat{H}}^{-1} L^\top )^{-1} L, \end{aligned}$$

(61)

where

$$\begin{aligned} L = \begin{bmatrix} e^{(u)}_E &{} -I_{v2u}^g e^{(v)}_W \\ c_1^2 I_{u2v}^b d_E^{(u)} &{} -c_2^2 d_W^{(v)} \end{bmatrix}. \end{aligned}$$

(62)

with the OP interpolation operators and

$$\begin{aligned} L = \begin{bmatrix} e^{(u)}_E &{} -I_{v2u} e^{(v)}_W \\ c_1^2 I_{u2v} d_E^{(u)} &{} -c_2^2 d_W^{(v)} \end{bmatrix}, \end{aligned}$$

(63)

with the traditional interpolation operators.

1.3 A3. Two-Block SBP-P-SAT

The two-block discretization with Neumann boundary conditions imposed using SBP-SAT and interface conditions imposed using the hybrid SBP-P-SAT method is given by

$$\begin{aligned} Q = P \left( \begin{bmatrix} c_1^2 D_L^{(u)} &{} 0 \\ 0 &{} c_2^2 D_L^{(v)} \end{bmatrix} + SAT_{IC} + SAT_{BC} \right) P, \end{aligned}$$

(64)

and

$$\begin{aligned} G(t) = -P \begin{bmatrix} {\bar{\tau }}_W^{(u)} {\textbf{g}}_W^{(u)}(t) + {\bar{\tau }}_S^{(u)} {\textbf{g}}_S^{(u)}(t) + {\bar{\tau }}_N^{(u)} {\textbf{g}}_N^{(u)}(t) &{} 0 \\ 0 &{} {\bar{\tau }}_E^{(v)} {\textbf{g}}_E^{(v)}(t) + {\bar{\tau }}_S^{(v)} {\textbf{g}}_S^{(v)}(t) + {\bar{\tau }}_N^{(v)} {\textbf{g}}_N^{(v)}(t) \end{bmatrix}. \end{aligned}$$

(65)

Here

$$\begin{aligned} SAT_{IC} = \begin{bmatrix} 0 &{} 0 \\ -(H_x^{(v)})^{-1} e^{(v) \top }_W c_1^2 I_{u2v}^b d_E^{(u)} &{} (H_x^{(v)})^{-1} e^{(v) \top }_W c_2^2 d_W^{(v)} \end{bmatrix}, \end{aligned}$$

(66)

with the OP interpolation operators and

$$\begin{aligned} SAT_{IC} = \begin{bmatrix} 0 &{} 0 \\ -(H_x^{(v)})^{-1} e^{(v) \top }_W c_1^2 I_{u2v} d_E^{(u)} &{} (H_x^{(v)})^{-1} e^{(v) \top }_W c_2^2 d_W^{(v)} \end{bmatrix}, \end{aligned}$$

(67)

with the traditional interpolation operators. Also

$$\begin{aligned} SAT_{BC} = \begin{bmatrix} {\bar{\tau }}_W^{(u)} d_W^{(u)} + \bar{\tau }_S^{(u)} d_S^{(u)} + {\bar{\tau }}_N^{(u)} d_N^{(u)} &{} 0 \\ 0 &{} \bar{\tau }_E^{(v)} d_E^{(v)} + {\bar{\tau }}_S^{(v)} d_S^{(v)} + \bar{\tau }_N^{(v)} d_N^{(v)} \end{bmatrix}, \end{aligned}$$

(68)

where the penalty matrices ${\bar{\tau }}^{(u)}_{W,S,N}$ and ${\bar{\tau }}^{(v)}_{E,S,N}$ are given by (57) evaluated in each respective block and ${\textbf{g}}^{(u)}_{W,S,N}$ and ${\textbf{g}}^{(v)}_{E,S,N}$ are boundary data. The projection operator is given by

$$\begin{aligned} P = I - {\hat{H}}^{-1} L^\top (L {\hat{H}}^{-1} L^\top )^{-1} L, \end{aligned}$$

(69)

where

$$\begin{aligned} L = \begin{bmatrix} e^{(u)}_E&-I_{v2u}^g e^{(v)}_W \end{bmatrix}. \end{aligned}$$

(70)

with the OP-interpolation operators and

$$\begin{aligned} L = \begin{bmatrix} e^{(u)}_E&-I_{v2u} e^{(v)}_W \end{bmatrix}, \end{aligned}$$

(71)

with the traditional interpolation operators.

1.4 A4. Two-Block SBP-SAT

The two-block discretization with Neumann boundary conditions imposed using SBP-SAT and interface conditions imposed using the SBP-SAT is given by

$$\begin{aligned} Q = \begin{bmatrix} c_1^2 D_L^{(u)} &{} 0 \\ 0 &{} c_2^2 D_L^{(v)} \end{bmatrix} + SAT_{IC} + SAT_{BC}, \end{aligned}$$

(72)

and

$$\begin{aligned} G(t) = -\begin{bmatrix} {\bar{\tau }}_W^{(u)} {\textbf{g}}_W^{(u)}(t) + {\bar{\tau }}_S^{(u)} {\textbf{g}}_S^{(u)}(t) + {\bar{\tau }}_N^{(u)} {\textbf{g}}_N^{(u)}(t) &{} 0 \\ 0 &{} {\bar{\tau }}_E^{(v)} {\textbf{g}}_E^{(v)}(t) + {\bar{\tau }}_S^{(v)} {\textbf{g}}_S^{(v)}(t) + {\bar{\tau }}_N^{(v)} {\textbf{g}}_N^{(v)}(t) \end{bmatrix}. \end{aligned}$$

(73)

The SAT imposing the interface conditions using the OP interpolation operators is given by

$$\begin{aligned} SAT_{IC} = \begin{bmatrix} H_x^{(u)} &{} 0 \\ 0 &{} H_x^{(v)} \end{bmatrix}^{-1} \begin{bmatrix} SAT_{IC}^{11} &{} SAT_{IC}^{12} \\ SAT_{IC}^{21} &{} SAT_{IC}^{22} \end{bmatrix}, \end{aligned}$$

(74)

where

$$\begin{aligned} \begin{aligned} SAT_{IC}^{11}&= SAT_{u1} e^{(u)}_E + SAT_{u2} I_{C2F}^g e^{(u)}_E + SAT_{\partial u} c_1^2 d_E^{(u)}, \\ SAT_{IC}^{12}&= -SAT_{u1} I_{F2C}^g e^{(v)}_W - SAT_{u2} e^{(v)}_W - SAT_{\partial u} c_2^2 I_{F2C}^b d_W^{(v)}, \\ SAT_{IC}^{21}&= -SAT_{v1} I_{C2F}^g e^{(u)}_E - SAT_{v2} e^{(u)}_E - SAT_{\partial v} c_1^2 I_{C2F}^b d_E^{(u)} \\ SAT_{IC}^{22}&= SAT_{v1} e^{(v)}_W + SAT_{v2} I_{F2C}^g e^{(v)}_W + SAT_{\partial v} c_2^2 d_W^{(v)} \end{aligned} \end{aligned}$$

(75)

and

$$\begin{aligned} \begin{aligned} SAT_{u1}&= -\frac{c_1^2 \tau }{h_x^{(u)}} e^{(u) \top }_E + \frac{c_1^2}{2} d_E^{(u) \top } \\ SAT_{u2}&= -\frac{c_2^2 \tau }{h_x^{(v)}} e^{(u) \top }_E I_{F2C}^b \\ SAT_{\partial u}&= - \frac{1}{2} e^{(u) \top }_E \\ SAT_{v1}&= -\frac{c_2^2 \tau }{h_x^{(v)}} e^{(v) \top }_W - \frac{c_2^2}{2} d_W^{(v) \top } \\ SAT_{v2}&= - \frac{c_1^2 \tau }{h_x^{(u)}} e^{(v) \top }_W I_{C2F}^b \\ SAT_{\partial v}&= \frac{1}{2} e^{(v) \top }_W \end{aligned}. \end{aligned}$$

(76)

Here $h_x^{(u,v)}$ are the grid sizes in the x-direction in each block,

$$\begin{aligned} \tau = \frac{\theta }{4 \alpha }, \end{aligned}$$

(77)

with $\theta = 3$ and $\alpha = 0.2508560249$ for the 4th order operators and $\alpha = 0.1878715026$ for the 6th order operators. With the traditional interpolation operators, we have

$$\begin{aligned} SAT_{IC} = \begin{bmatrix} H_x^{(u)} &{} 0 \\ 0 &{} H_x^{(v)} \end{bmatrix}^{-1} \begin{bmatrix} SAT_{IC}^{11} &{} SAT_{IC}^{12} \\ SAT_{IC}^{21} &{} SAT_{IC}^{22} \end{bmatrix}, \end{aligned}$$

(78)

where

$$\begin{aligned} \begin{aligned} SAT_{IC}^{11}&= SAT_{u1} e^{(u)}_E + SAT_{u2} I_{C2F} e^{(u)}_E + SAT_{\partial u} c_1^2 d_E^{(u)}, \\ SAT_{IC}^{12}&= -SAT_{u1} I_{F2C} e^{(v)}_W - SAT_{u2} e^{(v)}_W - SAT_{\partial u} c_2^2 I_{F2C} d_W^{(v)}, \\ SAT_{IC}^{21}&= -SAT_{v1} I_{C2F} e^{(u)}_E - SAT_{v2} e^{(u)}_E - SAT_{\partial v} c_1^2 I_{C2F} d_E^{(u)} \\ SAT_{IC}^{22}&= SAT_{v1} e^{(v)}_W + SAT_{v2} I_{F2C} e^{(v)}_W + SAT_{\partial v} c_2^2 d_W^{(v)} \end{aligned} \end{aligned}$$

(79)

and

$$\begin{aligned} \begin{aligned} SAT_{u1}&= - c_1^2 \tau e^{(u) \top }_E + \frac{c_1^2}{2} d_E^{(u) \top } \\ SAT_{u2}&= - \tau c_2^2 e^{(u) \top }_E I_{F2C} \\ SAT_{\partial u}&= - \frac{1}{2} e^{(u) \top }_E \\ SAT_{v1}&= - \tau c_2^2 e^{(v) \top }_W - \frac{c_2^2}{2} d_W^{(v)\top } \\ SAT_{v2}&= - \tau c_1^2 e^{(v) \top }_W I_{C2F} \\ SAT_{\partial v}&= \frac{1}{2} e^{(v) \top }_W \end{aligned}. \end{aligned}$$

(80)

Here

$$\begin{aligned} \tau = \theta \max \left( \frac{1}{\alpha h_x^{(u)}},\frac{1}{\alpha h_x^{(v)}} \right) \end{aligned}$$

(81)

with $\theta = 3$ and $\alpha = 0.2508560249$ for the 4th order operators and $\alpha = 0.1878715026$ for the 6th order operators. Here $h_x^{(u,v)}$ are the grid sizes in the x-direction in each block. The SAT imposing the boundary conditions is given by

$$\begin{aligned} SAT_{BC} = \begin{bmatrix} {\bar{\tau }}_W^{(u)} d_W^{(u)} + \bar{\tau }_S^{(u)} d_S^{(u)} + {\bar{\tau }}_N^{(u)} d_N^{(u)} &{} 0 \\ 0 &{} \bar{\tau }_E^{(v)} d_E^{(v)} + {\bar{\tau }}_S^{(v)} d_S^{(v)} + \bar{\tau }_N^{(v)} d_N^{(v)} \end{bmatrix}, \end{aligned}$$

(82)

where the penalty matrices ${\bar{\tau }}^{(u)}_{W,S,N}$ and ${\bar{\tau }}^{(v)}_{E,S,N}$ are given by (57) evaluated in each respective block and ${\textbf{g}}^{(u)}_{W,S,N}$ and ${\textbf{g}}^{(v)}_{E,S,N}$ are boundary data.