Boundary-optimized summation-by-parts operators for ﬁnite difference approximations of second derivatives with variable coeﬃcients

Boundary-optimized summation-by-parts (SBP) ﬁnite difference operators for second derivatives with variable coeﬃcients are presented. The operators achieve increased accuracy by utilizing non-equispaced grid points close to the boundaries of the grid. Using the optimized operators we formulate SBP schemes for the acoustic and elastic operators deﬁned directly on curvilinear multiblock domains. Numerical studies of the acoustic and elastic wave equations demonstrate that, compared to traditional SBP difference operators, the new operators provide increased accuracy for surface waves as well as block interfaces in multiblock grids. For instance, simulations of Rayleigh waves demonstrate that the boundary-optimized operators more than halve the runtime required for a given error tolerance. © 2023 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons .org /licenses /by /4 .0/).


Introduction
Since their initial development in the 1970s [17,18] and 1990s [36], high-order summation-by-parts (SBP) finite difference methods have emerged as an efficient and robust discretization method for hyperbolic partial differential equations.When paired with stable boundary treatment, such as simultaneous approximation terms (SAT) [6], projection [27], or ghost point methods [42], the SBP property facilitates the construction of energy or entropy stable schemes.For reviews of SBP finite difference methods paired with SAT we refer to [38,8].
Unless artificial dissipation is added to the scheme, SBP operators based on a diagonal norm are required to avoid stability issues for problems with variable coefficients, e.g., with variable material parameters or discretized on curvilinear grids [37,24].Moreover, schemes obtained utilizing such diagonal-norm difference operators naturally allow for explicit time integration methods and are thus well-suited for hyperbolic problems.However, the SBP property for diagonal-norm difference operators comes at the cost of a reduction in the order of accuracy in the boundary stencils.In particular, for a 2pth-order accurate diagonal SBP norm H, the associated SBP difference operator with interior order 2p is at most pthorder accurate at the boundary [20].To address the reduced boundary accuracy, boundary-optimized SBP operators for first derivatives were first presented in [25], and later improved upon in [26].These operators are also pth-order accurate near boundaries, but reduce the proportionality constant in the leading order error term significantly.This is achieved by utilizing stencils based on non-equispaced grid points close to the boundaries of the grid.For problems where boundary errors dominate, the reduction in error can be several orders of magnitudes compared to traditional equispaced SBP operators, as illustrated in [26].Moreover, while the traditional SBP difference operators of orders higher than 6 typically suffer from large truncation errors and therefore rarely are preferred in practice, the boundary-optimized operators of orders 8, 10 and 12 are shown to remain accurate.This allows for the use of coarser grids, and thus higher computational efficiency, to achieve a given error tolerance.
Currently, boundary-optimized SBP difference operators have been presented for first derivatives in [26], and for second derivatives with constant coefficients in [12].For second derivatives with variable coefficients, SBP operators with a minimal stencil width (so-called narrow-stencil operators) were presented in [23].These operators are defined on equispaced grids and exist for orders two, four and six.In this work, we present boundary-optimized narrow-stencil SBP difference operators for second derivatives with variable coefficients of even orders from 4 up to 12. Through numerical studies of the acoustic and elastic wave equations, we then demonstrate the benefit of using the new operators for problems on curvilinear multiblock domains, as well as problems involving surface waves.In particular, through a convergence and efficiency study of Rayleigh waves, we show that the novel boundary-optimized operators are preferable to the traditional narrow-stencil SBP operators when simulating surface waves in elastic media.
The article is developed as follows: In Section 2 we establish the notation used throughout the article.In Section 3 SBP operators are defined.Importantly, we outline the procedure used to construct the novel boundary-optimized variablecoefficient second derivative operators.In Section 4 we present discretizations of the acoustic and elastic operators on curvilinear multiblock domains, which we later make use of in the numerical studies presented in Sections 5 -6.Finally, the article is concluded in Section 7.

Notation
Throughout the article, boldface symbols will be used to denote discrete quantities, e.g. the discrete grid function approximating the function u(x), x ∈ is denoted u and is defined on the grid .Vector-valued functions and grid functions are denoted with a bar, e.g., the grid function restriction of ū is denoted ū.The approximation of scalar functions as diagonal operators is denoted with double bars, e.g. the operator approximating f (x) is f = diag(f).
The L 2 -inner product for scalar functions u, v on is denoted u, v = uvd . ( The notation used in the discrete setting is analogously for a quadrature rule H on .Similarly, define the bilinear form for integrating over a surface ∈ ∂ as u, v = uvd , (2.3) and in the discrete setting u, v = (e T u) T H (e T v), (2.4) where H is a quadrature rule on , and e are boundary restriction operators such that, e.g., e T ū ≈ u(x), x ∈ .
When formulating relations for vector-valued problems as well as when defining metric relations of coordinate transformations and operators on multiblock curvilinear grids, we will make use of index notation.Here Einstein's summation convention is applied to sum over repeated indices i, j, k, l and I, J , K , L used to denote spatial components.For instance, (2.1) and (2.3) may then be extended to vector-valued functions ū and v as In Section 4 we deal with transformations between a physical domain , and a reference domain ˜ .We denote the Cartesian coordinate axes of and ˜ by x, and ξ , respectively.Moreover, I, J , K , L will be used to denote indices corresponding to , whereas i, j, k, l refers to indices on ˜ .To denote a partial derivative on we write . The summation convention then also applies to partial derivatives.For instance the divergence of a vector-valued function ū is written as . Similarly, for a difference operator

SBP finite difference operators
We proceed by defining SBP finite difference operators in one space dimension, before extending them to multiple dimensions in Section 4. In particular, Section 3.3 presents the procedure used to construct the novel boundary-optimized second derivative operators.Consider therefore the discretization of the interval = [x l , x r ], by m grid points = [x 1 , . . ., x m ] T , where x 1 = x l , and x m = x r .The boundary-optimized SBP operators derived in [26] and the novel boundary-optimized second derivative operators presented in this article are defined for grids with non-equispaced grid points close to the boundaries, but equispaced interior grid points.Given n non-equispaced grid points at each boundary, and the interior grid spacing h, the n − 1 spacings in the left boundary region are h i = hd i , i = 1, . . ., n − 1, where d i is the grid spacing scaling for the subinterval [x i , x i+1 ].Moreover, the grid is symmetric about its center such that in the right boundary region d m−i = d i .See Fig. 3.1.For non-equispaced grids, the interior grid spacing is given by h = L 2D+(m−2n−1) , where D = n i=1 d i .The weights d i were derived in [26] and can be found therein, as well as in the online repository [35].On equispaced grids Finite difference grid with non-equispaced boundary grid points in the boundary regions.

Definitions
SBP operators are defined with respect to a discrete norm H, defining a quadrature rule on the grid.In this article, we exclusively consider diagonal H. Thus all the operators presented herein are so-called diagonal-norm SBP operators.In addition, we will make use of boundary restriction operators e l = [1, 0, . . .0] T , e r = [0, . . .0, 1] T .
The restriction operators simply restrict grid functions to the boundary grid point, e.g., e T l u = u(x 1 ).Note that although SBP discretizations may be generalized to grids in which the boundary points are excluded (so-called generalized SBP operators [7,9]), the operators presented herein include the boundary points.
First derivative SBP operators are given by the following definition [29]: ∂x is said to be a first derivative SBP operator, if for the norm H, the relation holds.
The focus in the present work is on variable-coefficient second derivative SBP operators compatible with Definition 3.1.Such operators of orders two, four and six, were first presented in [23,22]: ∂x , is said to be a compatible second derivative SBP operator, if for the norm H, the relation holds, where the first and last rows of D1 approximate ∂ ∂x and R(c) = R(c) T ≥ 0.
The remainder term R(c) consists of undivided difference approximations such that R(c) ≈ 0 up to the order of accuracy.
The structure of R(c) will be outlined in detail in Section 3.
In Section 4.3 we present the corresponding SBP properties for discretizations of the acoustic and elastic operators and briefly discuss how the relations may be used to construct energy-stable schemes.
In addition to Definition 3.2, the new boundary-optimized second derivative operators presented in this article, satisfy a stronger condition termed full compatibility [28].The following definition is therefore central to this work: Fully compatible D 2 (c) may therefore be obtained from the discretization D 1 cD 1 , together with an appropriate choice of the term R(c).This observation will be key when deriving D 2 (c), fully compatible with the boundary-optimized D 1 of [26].Another important observation is that Definition 3.3 holds iff the first and last row of D1 and D 1 in (3.2) are equal, such that D1 may be replaced by D 1 in Definition 3.2.Indeed, this is the definition of full compatibility used in [1,2].The difference between compatible and fully compatible operators is therefore in the boundary closure of D 2 (c).Importantly, full compatibility greatly simplifies the derivation of energy-stable schemes for wave equations, when Dirichlet-type conditions are imposed weakly [11].Lastly, we will consider second derivative operators with different stencil widths, and the following definition [29] is therefore of relevance: Definition 3.4.A finite difference operator using a 2pth-order accurate repeated interior stencil of minimal stencil width for the Cauchy problem is said to be a 2pth-order accurate narrow-stencil operator.
The operators of [29,23] are compatible diagonal-norm narrow-stencil operators defined for equispaced grids, and we will refer to these as the traditional SBP difference operators.Wide-stencil second derivative operators are constructed as which coincides with setting R(c) = 0 in Definition 3.2.Moreover, wide-stencil operators are also fully compatible, following Definition 3.3.

Accuracy and convergence rates
For equations with second derivatives in space, discretized with SBP operators with a 2pth-order accurate interior stencil and rth-order accurate boundary stencils, the expected spatial convergence rate of energy stable schemes is at least q = min(2p, r + 2) [14,43,39].As previously stated, diagonal-norm SBP operators with interior order 2p are at most pthorder accurate at the boundary.This holds for the traditional operators derived in [23], such that the expected asymptotic convergence rate is q = min(2p, p + 2).
For wide-stencil second derivatives, the boundary order is reduced to r = p − 1, and the expected asymptotic convergence rate then drops to at least q = min(2p, p + 1).As outlined in Section 3.3, the new boundary-optimized operators are adapted from a wide-stencil operator, such that the resulting operator is fully compatible.Such adapted fully compatible operators have boundary order r = p − 1.In practice, computations are most often performed on marginally resolved grids, and in such settings convergence rates above the asymptotic rate may be observed.In Sections 5 -6, we demonstrate that the novel boundary-optimized operators yield much smaller errors than the traditional operators on computationally tractable grids, despite the local reduction in the order of accuracy.To measure convergence rates we use the H-norm error, i.e., for e = u − u * , where u * is the grid function restriction of the true solution, the convergence rate is computed as where superscripts denote the number of grid points, m 1 , and m 2 , used in the computations.

Boundary-optimized SBP operators for second derivatives
In this section, the procedure to construct the novel boundary-optimized SBP operators for second derivatives with variable coefficients is presented.As previously mentioned, boundary-optimized SBP D 1 operators with non-uniform grid point distributions near boundaries were first presented in [25,26].The operators were derived by solving a non-linear optimization problem, where the parameters were the grid point locations and the coefficients of D 1 and H.The SBP property (Definition 3.1) was imposed as a constraint, along with the accuracy requirements that D 1 should be exact for polynomials of degree 2p in the interior and polynomials of degree p close to the boundaries.Any remaining free parameters were then optimized to minimize the two leading order error terms near boundaries (corresponding to polynomials of degree p + 1 and p + 2).Boundary-optimized SBP D 2 operators for constant-coefficient second derivatives of orders six and eight were derived in [12], following the same procedure as described below.We will now extend the operators to variable coefficients, and derive new boundary-optimized operators of orders 4, 10 and 12.
Consider a variable-coefficient second derivative SBP operator D 2 (c), fully compatible with the boundary-optimized D 1 of [26].Full compatibility implies that the D 1 and D 2 (c) operators must share the same grid and norm matrix H. Thus, we here consider both the grid point distribution and H as a priori determined by the D 1 operators.By Definition 3.3, a fully compatible second derivative operator is of the form D 2 (c) = D 1 cD 1 − H −1 R(c).As stated earlier, this means that a fully compatible operator may be obtained from a wide-stencil approximation D 1 cD 1 given an appropriate choice of the term H −1 R(c).Since D 1 is optimized for boundary accuracy, the boundary accuracy of the wide-stencil operator D 1 cD 1 is expected to be quite good (at least compared to the same operator obtained using traditional D 1 ).However, the wide-stencil approximation suffers from poor resolution characteristics for high-frequency modes, as illustrated in Section 3.4.Thus, the primary task at hand is to find R(c) such that D 2 (c) is a narrow-stencil operator.Further, R(c) should not degrade the boundary accuracy or cause a large spectral radius.
Given a 2pth-order accurate boundary-optimized operator D , we therefore start from the wide-stencil operator, , and construct R (2p) (c) such that the interior stencil of To obtain a minimal-width stencil we must cancel terms in the wide stencil while preserving the order of accuracy.This may be achieved by adding linear combinations of higher-order derivatives, scaled by the grid spacing to appropriate powers, such that (H From the above considerations, we make the ansatz: where the following constraints are to be met i is a narrow-stencil operator based on the grid point distribution of D , consistent with ∂ i ∂x i except close to the boundaries where it is allowed to be zero.i may be found in the online repository [35] where the new D 2 (c) operators, the boundary-optimized D 1 of [26] and their corresponding grids are available.

α (2p)
• Fourth order case (2p = 4): • Sixth order case (2p = 6): α (6)   1−3 = (80, 600, 3600) , ( • Eighth order case (2p = 8): α (8)   1−4 = (350, 2520, 14700, 78400) , ( C(8) Remark.The traditional variable-coefficient second derivative operators utilize pth-order accurate stencils in the boundary closure.These operators were constructed by requiring that Remark.It should be emphasized that the goal of this study is to derive compatible D 1 and D 2 (c) operators that can be combined to construct schemes for second-order PDEs with improved boundary accuracy.Here, we accomplish this by combining the existing boundary-optimized D 1 of [26] with a careful choice of R(c).However, there is no guarantee that the resulting combination of D 1 and D 2 (c) yields an optimal scheme for second-order PDEs.Fortunately, the numerical results presented in Sections 5 -6 demonstrate that the new operators do indeed provide higher accuracy than the traditional D 1 and D 2 (c) for a range of problems.An alternative, and perhaps more ambitious, approach would be to derive new compatible D 1 and D 2 (c) simultaneously, allowing for new grid-point distributions close to the boundaries.We foresee two major difficulties with this approach.Firstly, the resulting non-linear optimization problem is likely very challenging to solve.Secondly, the objective function in the optimization procedure would have to contain some combination of errors of D 1 and D 2 (c), and how to weight the different errors is not obvious.Nevertheless, this approach could be interesting to explore in a future study.

Resolution characteristics
In [12] boundary-optimized second derivative operators (with constant coefficients) were used in discretizing the incompressible Navier-Stokes equations.There it was observed that including two additional stencil points in the second derivative operators (compared to the minimal stencil width), resulted in a slightly more accurate scheme, while also allowing for a slightly larger time step in explicit time integration.Such semi-narrow-stencil operators may be obtained by omitting the lowest-order derivative term in (3.3).In the following, we illustrate that compared to wide-and semi-narrow-stencil operators, the narrow-stencil operators provide superior accuracy over the full length of scales that can be realized on the given grid.This property is termed the resolution characteristics of the approximation in [19] and is directly related to the dispersion error of the finite difference scheme for wave propagation problems.
To analyze the resolution characteristics we compute the symbols (i.e. the Fourier transformation) of the interior stencils of the SBP operators.The difference between the exact and discrete symbol then quantifies the resolution characteristics.This analysis (see for example [19,40]) is strictly speaking only valid for periodic domains, but it still says much about the interior stencil approximation on bounded domains.Consider a pure Fourier mode u = e ikx .Then ∂ 2 u ∂x 2 = −k 2 u.Consider the same Fourier mode on a grid over [−1, 1] with grid spacing h.It is convenient to introduce a scaled wavenumber κ = kh, where κ ∈ [0, π], such that we may write In the discrete setting, for a 2pth-order accurate centered finite difference stencil of width w, we correspondingly have D 2 u = D 2 u where the symbol D 2 can be written as and where s j are the stencil weights, with s 0 being the central element.See [15,19] for details.The difference P(κ) − κ 2 is then a measure of the error in the Fourier transformation, and defines the resolution characteristics of the difference approximation.In Fig. 3.2 we plot P(κ) for wide-, semi-narrow and narrow stencils of orders 2 to 12.As is well-known, the wide-stencil approximations have poor resolution characteristics for higher frequencies as shown in Fig. 3.2a.Of greater importance is the comparison of semi-narrow approximations and the corresponding narrow-stencil approximations presented in Fig. 3.2b.Clearly, the narrow-stencil approximations have superior resolution characteristics and are thus preferable for wave propagation problems.The semi-narrow stencil approximations are therefore not discussed further.

Wave operators on curvilinear multiblock grids
In this section, we introduce SBP discretizations of the acoustic operator (i.e., the variable-coefficient Laplace operator) and the elastic operator for curvilinear multiblock grids.SBP operators are extended to multi-dimensional curved domains using tensor products and coordinate transformations.For examples of this process applied to wave equations see, e.g., [33,30,1,2].By incorporating the metric transformation (in [3] referred to as 'encapsulation'), discrete operators with SBP properties may be defined directly on the curvilinear domain.Similarly, one may also incorporate the interfaces of multiple connected curvilinear domains (see, e.g., [2,21]), and in this way construct multiblock SBP operators on the entire physical domain.
Throughout the article, we restrict ourselves to 2D problems, but extensions to higher dimensions follow straightforwardly.See in particular [1,2] where the derivation of the acoustic and elastic operators on curvilinear multiblock grids are presented in detail for multiple dimensions.In this article, the analysis is restricted to fully compatible operators, to simplify boundary and interface treatment.As mentioned in Section 2, we will make heavy use of Einstein's summation convention.The presentation follows that of [1,2] fairly closely.Moreover, the following metric relation holds [41] J

Coordinate transformation
Next, we relate the normals n and ν.To do this, we first introduce the surface scale factor J , which in 2D is defined as where Next, let contain an elastic medium with stiffness tensor C I J K L .Then the elastic operator is given by ∂ I C I J K L ∂ K .Following [2], we define the transformed stiffness tensor on ˜ as Cijkl = F Ii J C I jKl F Kk .Using (4.1) and ( 4.2), one can show that the elastic operator satisfies Moreover, we introduce the traction operator T J L as such that the traction vector τ is given by τ J = T J L u L .
In the following section, we will use equations (4.5), and (4.7) to construct SBP approximations of the wave operators based on difference operators defined on the reference domain.

Curvilinear operators
To define SBP operators for the non-rectangular grid , we start by defining the operators on ˜ and then make use of the above metric relations.Consider therefore the grid ˜ = ξ 1 ⊗ ξ 2 discretizing ˜ = [0, 1] 2 using m ξ 1 × m ξ 2 grid points.We use the notation A ξ i to denote a 1D operator A on the grid ξ i .Operators not dependent on variable coefficients are conveniently extended to multiple dimensions using tensor products.For non-variable operators A ξ i and identity operators From here on, we will simply use the subscript i in place of ξ i , e.g., D i := D ξ i , and similar for subscripts j, k, l.
Constructing a variable-coefficient narrow-stencil second derivative operator on ˜ requires slightly more care.An opera- tor for the mixed derivative ∂ i c(x)∂ j , i = j, is straightforwardly constructed as D i cD j .However for i = j, the operator D ii (c) (no sum over i) is given by applying the operator D ξ i 2 (c) of Definition 3.2 for each grid line.Thus we write (4.10) Next, the metric coefficients J , F Ii and J are approximated using difference approximations.To start with, note that F Ii can be expressed in terms of FiI , since F = F −1 by the inverse function theorem.In 2D the resulting relations are Thus by setting FiI = D i x I , we have J = F11 F22 − F12 F21 , and (4.11) is discretized using J and FiI to obtain FIi .Using J and FIi the remaining discrete metric quantities are straightforwardly obtained.For instance, a discrete first derivative on is given by discretizing (4.1), such that where D ij is given by (4.10) and α ij is the approximation of (4.6).Expanding the summation convention for 2D and using (4.10), (4.13) reads Similarly, the elastic operator ) is given by where again D ik is given by (4.10), and Ci JkL = FIi J CI J K L F Kk is a discretization of the continuous transformed stiffness tensor Ci JkL .
We conclude the section by defining discrete inner products and bilinear forms on .An SBP norm for is given by H = JH ˜ .
Note that since J approximates the volume scale factor J , H is a quadrature rule for such that indeed u, v ≈ u, v .
The SBP quadrature for a boundary is given by where J is the discretization of the surface scale factor J in (4.3), and H ˜ is the SBP quadrature for ˜ .For instance, on the boundary ˜ W , H ˜ = H ξ 2 .A bilinear form for the entire boundary ∂ is then given by

Summation-by-parts properties
Integration by parts is key in deriving energy balances and showing conservation properties of the continuous equations, and SBP operators allow us to mimic the analysis in the discrete setting.It is shown in [1,2] that the discrete wave operators of (4.13) and (4.15) satisfy SBP properties on , consistent with the continuous counterparts.The results are reiterated here for completeness, although proofs are omitted.We start by stating the continuous integration-by-parts properties, before stating the corresponding SBP property.
For scalar functions u and v, the acoustic operator (κ) satisfies the integration-by-parts property v, (κ)u = v, κ∂ nu ∂ − ∂ I u, κ∂ I v , (4.16) where ∂ n = nI ∂ I .In the discrete setting, we start by defining the discrete normal derivative as D n = nI D I where n is obtained by discretizing (4.4), and D I is given by (4.12).Then, the discrete acoustic operator in (4.13) satisfy the SBP property (4.17)where R = d i=1 R i ( κα ii ).Note that R = R T ≥ 0, and that R ≈ 0. Thus (4.17) is a discrete counterpart to the integrationby-parts property (4.16).
Similarly, for vector-valued functions ū, v, the elastic operator ∂ I C I J K L ∂ K satisfies the integration-by-parts property

.18)
A discrete traction operator is given by discretizing (4.8), resulting in Then, the discrete elastic operator in (4.15) satisfies the SBP property where W J L = k R k ( Ck JkL ).Due to symmetry properties of the stiffness tensor and R k , it follows that W J L = W T J L ≥ 0.Moreover, since R k ≈ 0 to the order of accuracy, it follows that W J L ≈ 0. Thus (4.19) mimicks the integration-by-parts relation (4.18).
The SBP properties (4.17), and (4.19) allow for the derivation of discrete energy rates for the acoustic and elastic wave equations respectively.This is achieved by replacing v and v by u and u, in which case the volume terms •, • , contribute to the discrete potential energy rate.The discrete energy rates are then useful for guiding the enforcement of boundary conditions, such that an energy stable scheme may be obtained.See, e.g., [34,10,1,2].

Multiblock operators
For problems discretized on multiblock grids, it is convenient to incorporate the interface treatment into the discrete wave operators.In this way, an SBP operator is defined for the entire multiblock grid.We will consider the special case of two connected domains.The generalization to multiple connected domains follows in the same fashion.Consider therefore the curvilinear multiblock domain = − ∪ + in Fig. 4.2, and denote the interface by = − ∩ + .We will use superscripts ± to indicate quantities defined on the respective domains.The two curvilinear domains are discretized such that ± conform along .In addition, we assume that the surface scale factors J ± coincide on the interface such that grid points on the reference grids ˜ ± conform along ˜ .We start by considering the acoustic operator.On each domain we discretize ± (κ ± )u ± as D ± (κ ± )u ± according to (4.13).Moreover, we impose the following interface conditions Note that it is not required that κ − = κ + on .Interface penalties SAT − and SAT + , enforcing (4.20) on − and + respectively, are given by SAT where with H 11 being the first element in the 1D SBP norms H ξ i .The parameter γ I specifies the strength of the penalty term.Broadly speaking, increasing γ I leads to stronger enforcement of the interface condition, at the cost of additional stiffness.Typically one chooses γ I close to the stability limit.Note that (4.22) is identical to the borrowing term given in [1] for fully compatible operators.
Define grid functions on the multiblock domain Then the discretization of (κ)u, on with interface conditions (4.20) is This leads to the formulation of the multiblock acoustic operator as (4.23)where Next, we formulate the SBP property for D M B (κ).
. Moreover, we define the interface jump operator • as Aw = A + w + − A − w − for A ± , w ± on .Then, for grid functions u, v, approximating continuous u, v on , satisfying (4.20), the SBP property of where The property is obtained by applying (4.17) to D ± (κ ± ).Note that the term Z A (u, v) ≈ 0 to the order of accuracy of the discretization, since R M B ≈ 0, and u ≈ v ≈ 0. Thus (4.24) is a consistent and accurate approximation of (4.16).For details see [1].
The elastic multiblock operator is formulated similarly.In the elastic setting, the interface conditions are The discrete elastic multiblock operator with interface conditions enforced using SAT is then given by where now the interface operators S (1)± J L , S (2)± J L are given by and (4.27) By using (4.19) one can show that the elastic multiblock operator satisfies the following SBP property (see [2]): where Again, the remainder term Z E ( ū, v) ≈ 0 to the order of accuracy of the scheme, for any smooth grid functions ū, v consistent with the interface conditions (4.25  [1] and Theorem 6 of [2].Moreover, one can show that the operators are self-adjoint.Self-adjointness is derived by applying the SBP property (4.17) to the right-hand-side of (4.24) once more.See Theorem 7 of [2] for a proof of self-adjointness of the elastic operator.The proof of self-adjointness of the acoustic operator follows analogously.Throughout the article, we use γ I = 1, which is right on the limit of stability.

Waves on a disk
To illustrate the benefit of using boundary-optimized SBP operators for curvilinear multiblock domains, consider the acoustic wave equation in a disk = {x = (x 1 , x 2 ) : |x| ≤ 1} with a Dirac delta source placed in the origin.Given zero initial conditions and a Dirichlet boundary condition u = g, the problem is thus (5.1) The time-dependency f (t) of the source is chosen to be a Gaussian centered about the time t s , i.e., (5. 2) The free space solution to (5.1) is obtained by changing to polar coordinates, applying Fourier transforms in time and Hankel transforms in space [5] and is given by dω, (5.3)where r = |x|.The integrand of (5.3) vanishes quickly as ω increases and thus the integral can be computed numerically.This is done using the Matlab function integral, specifying a tolerance of 10 −10 .Thus we solve (5.1) using (5.3) to set the Dirichlet boundary data g(x, t), up to some final time t = T before the wave has exited the domain.
The domain is discretized using five curvilinear blocks, as illustrated in Fig. 5.1a, where the grid point distribution of the eighth-order accurate boundary-optimized SBP operators is used.Discretizing in space using (4.23) (generalized to five blocks) results in where δ is a discretization of the Dirac delta function, and the SAT imposes the Dirichlet conditions in (5.1).Let b denote one of the four boundary segments of the outer blocks of .For the boundary-optimized fully compatible operators, the SAT is then where g b is the grid function restriction of the Dirichlet boundary data to b .The condition γ D ≥ 1 on the penalty parameter is required for energy stability.We use γ D = 1, which is on the limit of stability.For the traditional operators, we use the non-stiff SAT presented in [1].The Dirac delta function δ is discretized as in [31] using the same number of moment conditions as the order of accuracy of the SBP operators.For wide-stencil operators, additional smoothness conditions are required on the source discretization for the scheme to be accurate.In the wide-stencil computations, the same number of smoothness and moment conditions are therefore used.
The problem is solved with wave speed c = 1 up to t = 1 using source parameters σ = 1/25, and t s = 0.3, by timeintegrating (5.4) using the standard fourth-order Runge-Kutta method.At the final time t = 1, the wave has interacted with all the interfaces of the multiblock grid.See Fig. 5.1b.The time step is chosen as k = 0.01h (where h is the minimal grid  spacing on the multiblock grid), which is small enough for spatial errors to dominate.Errors and convergence rates using traditional and boundary-optimized narrow-stencil operators are presented in Tables 5.1-5.2where m denotes the number of grid points along the x 1 and x 2 coordinate axes.Comparing the fourth-, and sixth-order accurate operators, it is clear that the boundary-optimized operators yield smaller errors.For the fourth-order operators the gain in accuracy is slight, while for sixth-order operators, the gain in accuracy is approximately 1 magnitude on the coarser grids.The 8th-order operators perform particularly well, while for orders 10 and 12, the errors on coarser grids are relatively large and the observed high convergence rates are pre-asymptotic.Table 5.3 presents the results obtained using wide-stencil boundary-optimized operators.In addition to the narrow-stencil operators being computationally more efficient due to the smaller interior stencil widths, the improved dispersive properties of the narrow-stencil operators result in lower errors for all orders.In particular, for the fourth-and sixth-order operators, the boundary-optimized wide-stencil operators are less accurate than the corresponding traditional narrow-stencil operators.

Spectrum comparison
In this section we study the eigenvalue spectra of spatial discretizations of the wave equation, comparing the traditional and boundary-optimized operators for Dirichlet, Neumann, and interface conditions.Of particular interest is the spectral radius of the spatial discretization, since it affects the maximal stable time step allowed for explicit time integrators.For simplicity, the study is restricted to a 1D Cartesian grid on the domain = [0, 1] with the wave speed c = 1.When evaluating the interface coupling we impose the interface conditions (4.20) at each of the boundaries, such that the left and right boundaries are connected.The spatial discretization is then given by D = D 2 + SAT, where SAT differs depending on the boundary condition imposed.The SAT for Dirichlet conditions is the 1D version of (5.5) for the boundary-optimized operators.Similarly, for interface conditions, the boundary-optimized operators use the 1D version of (4.21).For the traditional operators, the Dirichlet and interface SATs can be found in [2].For Neumann conditions, we simply cancel the boundary terms in (3.2) To remove the dependency on the grid resolution, we consider the eigenvalue spectrum λ(h 2 D).The corresponding spectral radius is then given by ρ( ) is presented in Tables 5.4 -5.5.The results show that the boundary-optimized operators have larger spectral radii, in particular when Neumann conditions are imposed, in which case the spectral radii are approximately twice as large for comparable orders of accuracy.For the interface condition, the spectral radius is dictated by the largest spectral radii of  the Dirichlet and Neumann cases.For explicit time integrators, the CFL condition dictates that the time step is chosen as Thus the traditional operators should allow for larger time steps.For instance, comparing ρ(h 2 D) for the sixth-order accurate case with Neumann conditions, the traditional operators allow for approximately a 1.5 times larger time step.Although this is a significant restriction of the time step for the boundary-optimized operators, the numerical studies in Section 6 show that the increased boundary accuracy still pays off in terms of efficiency, in that coarser grids may be used to achieve a specified error tolerance, resulting in shorter runtimes.The eigenvalue spectra of the fourth-and sixth-order operators are presented in Fig. 5.2.

The elastic wave equation
Problems involving surface waves are another setting in which increased boundary accuracy is of interest.To illustrate the benefits of the new boundary-optimized operators we will therefore solve the elastic wave equation in settings where surface waves are present.The section begins by focusing purely on surface wave modes by simulating Rayleigh waves.In Section 6.2 elastic waves in a piecewise uniform isotropic medium excited by a surface source are studied.

Consider the elastic wave equation in the domain
where u J is the displacement in the x J -direction.A free surface condition is imposed on the northern boundary N .On the west and east boundaries, W and E , periodic conditions are imposed.The initial conditions and the traction boundary condition on the south boundary S will soon be discussed.We consider an isotropic material, where density and the first and second Lamé parameters are all unity, i.e., ρ = 1, λ = 1 and μ = 1.In this setting, the stiffness tensor is given by This results in pressure wave speed c P = √ 3, and shear wave speed c S = 1.For an open half-space in the x 2 -direction, the problem admits a solution of the following form: where the Rayleigh wave speed is c R = 2 − 2 √ 3 c S , A is a constant, χ is the wave number and See Chapter 6.1 of [13] for details.As evident from (6.3) and illustrated in Fig. 6.1, u 1 and u 2 are of maximal magnitudes at the surface x 2 = 0 and decrease exponentially as x 2 decreases.Thus the majority of the particle motion is located close to the surface boundary N and the problem should therefore be a good fit for testing out the new boundary-optimized operators.Moreover, the solution is periodic in the x 1 -direction, such that the only required boundaries are N and S .Discretizing (6.1) in space using (4.15) results in: (6.4) The SAT imposes the traction conditions on N and S , and is given by (see [10,2]): SAT J = SAT N J + SAT S J , (6.5) SAT S J = −H −1 e S H S (e T S T J L u L − g J ). (6.7) The initial conditions u 0 and u0 , and the boundary data g J (x, t) at the south boundary S are specified using (6.3).In the x 1 -direction periodic operators are used (i.e., the narrow interior stencil is applied for all grid points), and no boundary treatment is needed.In the x 2 -direction we use narrow-stencil SBP operators.
The problem is solved to t = 1, setting L 1 = L 2 = 4, A = 0.01 and χ = 4π .The domain is discretized using (m − 1) × m grid points.Time integration of (6.4) is performed using a sixth-order accurate explicit Runge-Kutta method, with the time step k = 1 10c P h.We compare the error obtained using the boundary-optimized operators, and the traditional adapted fully compatible operators presented in [2].The latter operators have reduced accuracy in a single grid point on each boundary, and an asymptotic convergence rate of q = p + 1.5 was observed in [2].Tables 6.1 -6.2 present the obtained errors and convergence rates.Similar to the results in Section 5, the results show that the boundary-optimized operators provide increased accuracy for a given grid resolution.
Next, we illustrate the efficiency gain in using boundary-optimized operators.This is done by comparing the error against the runtime while using a time step k close to the maximally stable time step.Consider therefore the same setup, but now using fourth-order Runge-Kutta method, where we choose the time steps as  imposed on the boundaries.These are Neumann-type conditions, for which the spectral radius is smaller for the traditional operators compared to the boundary-optimized operators (as shown in Section 5.2).Therefore the traditional operators allow for larger timesteps.The results of Tables 6.3 -6.4 align well with the obtained spectral radii in Tables 5.4 -5.5.In Fig. 6.2 the error is plotted against the runtime, taken as the minimum over 10 runs, using m = [31,36,41,46,51,61,71,81,91].
The computations are performed on a 2,6 GHz 6-Core Intel Core i7 processor, with 16 GB 2400 MHz DDR4 memory, using Matlab 2021a.The results show that above an error threshold of e = 5 • 10 −3 , the fourth-and sixth-order operators are fairly comparable in terms of efficiency, with the fourth-order boundary-optimized operator performing slightly better.Below this threshold the boundary-optimized operators provide better efficiency, with the 8th-order operators performing particularly well, until e = 6 • 10 −5 where the 10th-order operator becomes more beneficial to use.For instance, at the error threshold e = 1 • 10 −3 , the runtime using the sixth-order accurate traditional operators is approximately 0.2625 s, while the sixth-and eighth-order accurate boundary-optimized operators have a runtime of approximately 0.1703 s and 0.1448 s, respectively.The speedup obtained changing from traditional to boundary-optimized operators is therefore a factor  are stated in Tables 6.5 -6.6.The efficiency results are presented in Fig. 6.3.This time the results show that the boundaryoptimized operators always are beneficial to use.For instance, on the coarsest grid the runtimes of the fourth-order accurate operators are comparable.However, to achieve the same error as the fourth-order boundary-optimized operators using traditional operators, the runtime increase from increased resolution is approximately a factor of 3.33, i.e., the speedup changing from traditional to boundary-optimized operators is in this case 3.33.To achieve the error e = 4.6 • 10 −3 , i.e., the error on the finest grid using 6th-order traditional operators, changing to 10th-order boundary-optimized operators results in a speedup of approximately 6.13.

Excitations from a surface source in a piecewise uniform isotropic medium
Consider a half-space consisting of two vertically stacked elastic isotropic media, where the domain of interest is given by = + ∪ − , illustrated in Fig. 6.4.The density is constant throughout , but the first and second Lamé parameters differ in + and − , such that there is a jump in the stiffness tensor across the interface = + ∩ − .The material parameters are ρ + = 2.65, λ + = 0.0205, μ + = 2.24, ρ − = 2.65, λ − = 30.57,μ + = 20, and the stiffness tensor for each block is computed according to (6.2).The medium is initially at rest.At the surface N (i.e., the north boundary of + ), a free surface condition is imposed.The other boundaries of are artificial.Moreover, a time- together with the interface conditions (4.25) on .The time-dependency of the source term is a Ricker-wavelet centered about time t s , with peak frequency α P , and is given by where A is a constant.In the numerical experiments, we will use A = 0.1, α P = 10 and t s = α −1 P .Moreover, α P is used to define a measure of points per wavelength (PPWL).This is done as in [2], where the highest frequency that needs to be resolved on the grid is defined as the frequency α M for which the amplitude spectrum of f (t) is 5% of peak amplitude.As discussed in [2], the definition yields α M ≈ 2.4α P .We then define PPWL = min (c S )

2.4α P h 1
, where h 1 is the grid spacing in the x 1 -direction.
The problem is discretized in space using (4.26) with the free surface condition imposed using the SAT (6.6).To truncate the domain of interest, super-grid scale absorbing layers are used.The artificial boundaries of ± are thus connected to super-grid regions through numerical interfaces.In the super-grid regions, waves are damped out by a combination of  grid-stretching and artificial dissipation.The technique is outlined in [4,32,16,2] and is not addressed here further.The semi-discrete problem in the domain of interest is thus where D M B I K (C I J K L ) also incorporates the super-grid interfaces.The Dirac delta source is discretized as in [31], with the same number of moment conditions as the order of accuracy of the SBP operators, and zero smoothness conditions.Using m × ((m − 1)/2 + 1) grid points in each of ± , and integrating (6.8) using the fourth-order Runge-Kutta method with the time-step k = 1 10 max (c P ) h, the problem is solved up to t = 2.A reference solution is constructed using the sixth-order accurate traditional operators, and m = 321, corresponding to PPWL = 12.26.We then compute the surface displacement at the receiver on coarser grids with m = 81, m = 121, and m = 161, corresponding to PPWL = 3.06, PPWL = 4.60 and PPWL = 6.13 respectively.The results using the best-performing operators, namely the eighth-order accurate boundary-optimized operators, and sixth-order accurate traditional operators, are presented in Fig. 6.5.Overall, the boundary-optimized operators are more accurate in capturing the initial peaks between t = 0.5 and t = 1, corresponding to the arrival of surface waves.Moreover, on the coarser grids, u 2 obtained with traditional operators exhibits oscillations around t = 1.25 which are not present when using the boundary-optimized operators.At PPWL = 3.06, both solutions are fairly underresolved, while at PPWL = 6.13 the solution using boundary-optimized operators fits the reference solution very well.Displacement in the entire domain at times t = 0.25, 0.66, 1.32 are shown in Fig. 6.6.The solution is obtained using the eighth-order accurate boundary-optimized operators with PPWL = 6.13.At t = 0.25 waves are seen propagating from the source.At t = 0.66 surface waves have just reached the receiver, and reflections against and transmissions into the lower part of the domain are visible.At t = 1.37 the surface waves have exited the domain, and waves remaining in the domain are mainly due to the interaction between the interface and the free surface.
Remark.For this problem, the boundary-optimized and traditional operators of orders four and six performed similarly.For the boundary-optimized operators of orders 10 and 12, the solutions were underresolved on the coarser grids, similar to what was observed in Section 5.These results are omitted for the sake of brevity.

Conclusion
Boundary-optimized SBP operators for second derivatives with variable coefficients of orders 4, 6, 8, 10, and 12 have been developed.The performance of the operators in terms of accuracy and efficiency has been evaluated by solving the acoustic and elastic wave equation.The new operators have a minimal interior stencil width and are fully compatible with the boundary-optimized first derivative operators of [26].The full compatibility condition simplifies weak enforcement of Dirichlet-type conditions for second-order wave equations.The operators are available online at [35].
The operators provide higher accuracy for problems where boundary or interface errors are dominating when compared to the traditional operators of [23].The increased accuracy allows for the use of coarser grids to achieve a given error tolerance, resulting in lower computational times.In our studies, the use of eighth-order accurate boundary-optimized operators seems to give a balanced tradeoff between accuracy and computational cost.
In a future study, we will develop operators for coupling boundary-optimized finite difference schemes across nonconforming grid interfaces.Such coupling operators are non-trivial to design since they must be both highly accurate and handle the irregular grid point distribution of the non-conforming grid interfaces arising from the boundary-optimized operators.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
u I , v I = d I=1 u I v I d , u I , v I = d I=1 u I v I dwhere d is the dimension of .The summation convention also applies to discrete relations, such that for grid functions v and ū the extension of (2.2) and (2.4) is u I , v I = d I=1 u T I Hv I , u I , v I = d I=1 (e T u I ) T H (e T v I ).

Definition 3 . 3 .
Let D 1 and D 2 (c) be compatible first and second derivative SBP operators.If in addition, D 2 (c) = D 1 cD 1 − H −1 R(c) holds, the operators are said to be fully compatible.

2 () 2 (
c) is of minimal width.By construction R (2p) (c) is symmetric, and Constraint 1 further ensures semi-definiteness.Moreover, it implies that C(2p) i is an identity matrix for c(x) = 1, in which case R 2p (c) coincides with the constant coefficient case presented in [12].Constraints 2 -3 are both required for D (2pc) to obtain a 2pth-order accurate minimal-width interior stencil.Constraint 2 further ensures that R (2p) (c)u ≈ 0 to at least order p at the boundary.This can be deduced from the consistency of D (2p) p+ j , together with the fact that the terms h 2p−1+2 j (D (2p) p+ j ) T are O(h p ).The boundary accuracy of D (2p) 2 is then at least p − 1, which is consistent with that of the wide-stencil operator.Moreover, by allowing for D (2p) i to be zero close to the boundary, the spectral radius of D (2p) 2 (c) may be reduced.Through (3.3) and the above constraints we have constructed narrow-stencil D (2p) 2 (c) for orders 2p = 4, 6, 8, 10, 12. Below we present α (2p) j and C(2p) i at the interior points, for 2p = 4, 6, 8.The cases 2p = 10, 12 can be found in the appendix.At the boundaries, we use the same weights to compute the diagonal entries in C(2p) i , but with skewed stencils.The exact form of at the boundaries and then choosing C(2p) i in a careful way such that the boundary accuracy is preserved.It is currently not clear if one may find such a combination of D boundary-optimized operators, or indeed if boundary-optimized second derivative operators with pth-order accurate boundary closures even exists for the given grid point distribution.
Consider a smooth one-to-one mapping (ξ 1 , ξ 2 ) ∈ ˜ → (x 1 , x 2 ) ∈ between the physical domain and a square reference domain ˜ = [0, 1] 2 , illustrated in Fig. 4.1.Denote the boundaries by ∂ = { S , E , N , W }, and ∂ ˜ = { ˜ S , ˜ E , ˜ N , ˜ W }. Let F denote the Jacobian of the transformation, i.e., FiI = ∂x I ∂ξ i , and let J = det( F ) > 0 be the volume scale factor of the transformation.Similarly, let F be the Jacobian of the inverse transformation, i.e., F Ii = ∂ξ i ∂x I .As stated in Section 2 the notation ∂ I ≈ ∂ ∂x I and ∂ i ≈ ∂ ∂ξ i is used together with the summation convention.By the chain rule, the first derivative satisfies∂ I = F Ii ∂ i .

(4. 3 )
This allows us to relate the normals n and ν by Nanson's formula nI = J −1 J F Ii νi .

(4. 4 )
Now we are ready to present the transformation between and ˜ for the wave operators considered in this study.Let contain an acoustic medium with bulk modulus κ.Then the acoustic operator is simply the variable coefficient Laplace operator (κ) = ∂ I κ∂ I .By (4.1) and (4.2) the variable coefficient Laplace operator satisfies

(4. 9 )
Extensions to higher dimensions follow analogously.In this way the operators H ξ i and D ξ i on ˜ are constructed by extending the 1D operators H ξ i and D ξ i 1 from Section 3. Similarly, boundary restriction operators e ˜ , where ˜ denotes one of the boundaries of ˜ , are constructed by extending e ξ i l , and e ξ i r .Moreover, a norm for ˜ is given by

( 4 . 12 )
Now we are ready to define the discrete acoustic and elastic operators on .The acoustic operator D (κ) ≈ (κ) in (4.5) on is given by

Fig. 5 . 1 .
Fig. 5.1.(a) Multiblock grid using the grid point distribution of the eighth-order accurate boundary-optimized SBP operators.(b) Analytic solution (5.3) at t = 1.The wave has interacted with all of the multiblock interfaces.(For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

m f − 1 m− 1 k
(2p) f h.Here m f is the finest grid size, and k (2p) f is the corresponding largest stable time step obtained experimentally using 2pth-order accurate operators.We are particularly concerned with the accuracy on coarser grids and therefore focus on the range m ∈ [31, 91], i.e., m f = 91.In Tables 6.3 -6.4,k (2p) f is listed for the traditional and boundary-optimized operators.Note that traction conditions are
).Assuming proper enforcement of boundary conditions, the acoustic and elastic multiblock operators (4.23) and (4.26) are energy stable provided γ I ≥ 1 in (4.22) and (4.27).As discussed at the end of Section 4.3, energy stability is derived by setting v = u in (4.24) and v = u in (4.28), in which case Z A (u, u), and Z E ( ū, u) can be included in the respective potential energy rates.See for instance Theorem 3 of

Table 5 .1 Errors
e and rates q for (5.4) at t = 1, using traditional SBP operators.m l o g 10 ( e ) q l o g 10 ( e )

Table 5 .
q l o g 10 ( e ) q l o g 10 ( e ) q l o g 10 ( e )

Table 6 . 1
Errors e and rates q for the isotropic Rayleigh wave problem at t = 1, using traditional adapted SBP operators.

Table 6 . 3
Maximal stable time steps obtained for m f = 91 and L 1 = L 2 = 4, using traditional operators.

Table 6 . 6
Maximal stable time steps obtained for m f = 91 andL 1 = 4, L 2 = 8, usingTo provoke boundary errors further, the domain in the x 2 -direction is increased, setting L 2 = 8, while keeping the number of grid points fixed.Thus the resolution in x 2 is decreased.Moreover, the maximal stable time step is limited by the stencil in the x 1 -direction, and k f is similar for the traditional and boundary-optimized operators.The obtained values for k