arXiv:1311.4984v1 [math.NA] 20 Nov 2013

REVIEW OF SUMMATION-BY-PARTS SCHEMES FOR INITIAL-BOUNDARY-VALUE PROBLEMS
MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M
Abstract. High-order ﬁnite diﬀerence methods are eﬃcient, easy to program, scales well in multiple dimensions and can be modiﬁed locally for various reasons (such as shock treatment for example). The main drawback have been the complicated and sometimes even mysterious stability treatment at boundaries and interfaces required for a stable scheme. The research on summation-byparts operators and weak boundary conditions during the last 20 years have removed this drawback and now reached a mature state. It is now possible to construct stable and high order accurate multi-block ﬁnite diﬀerence schemes in a systematic building-block-like manner. In this paper we will review this development, point out the main contributions and speculate about the next lines of research in this area.
Keywords: well posed problems, energy estimates, ﬁnite diﬀerence; ﬁnite volume; boundary conditions; interface conditions; stability; high order of accuracy
1. Introduction
The research on Summation-By-Parts (SBP) schemes was originally driven by applications in ﬂow problems, including turbulence and wave propagation. The objective was to use highly accurate schemes to allow waves and other small features to travel long distances, or persist for long times. One of the ground-breaking papers showing the beneﬁt of high-order ﬁnite-diﬀerence methods for wave propagation problems is [KO72]. However, it has until recently proven diﬃcult to show the same beneﬁt in realistic simulations. Although it is easy to derive high-order ﬁnite diﬀerence schemes in the interior of the domain it is non-trivial to ﬁnd accurate and stable schemes close to boundaries. Furthermore, complicated geometries necessitates multi-block techniques. This poses yet another challenge for high-order ﬁnite diﬀerence schemes since solutions in diﬀerent blocks must be glued together in a stable and accurate way. The stencils near boundaries and block interfaces create diﬃculties. We will focus on the so called Simultaneous-Approximation-Term (SAT) technique where the boundary and interface conditions are imposed weakly.
The fundamental idea of SBP-SAT schemes is to allow proofs of convergence for linear and linearized problems. Convergence proofs form the bedrock of numerical analysis of PDEs since they provide the mathematical foundation that gives credibility to a numerical simulation. Without a proof of convergence, there is no guarantee that the numerical solution has any value at all. The conﬁdence that a discrete solution is an approximation of the true matematical solution is crucial, not only in practical engineering simulations, but also for the possibility to evaluate the accuracy of the model (i.e. the governing equation) itself and propose improved models. Without this quality assurance it is impossible to distinguish between modeling errors and numerical errors.
The SBP ﬁnite diﬀerence operators were ﬁrst derived in [KS74, KS77] and approximate coeﬃcients calculated. In [Str94], the analysis was revisited and exact
Date: November 21, 2013. 1

2

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

expressions for the ﬁnite diﬀerence coeﬃcients were obtained. However, the SBP ﬁnite diﬀerence operators alone, only admits stability proofs for very simple problems, and the use was limited. This changed with [CGA94] when Simultaneous Approximation Terms were proposed to augment SBP schemes. These are penaltylike terms that enforces boundary conditions weakly. With both SBP operators and the SAT technique at hand, stability proofs for more complicated systems of partial diﬀerential equations (PDEs) were within reach.
Finite diﬀerence methods are by no means the only choice of high-order schemes. There are numerous other high-order methods with diﬀerent strengths and weaknesses. However, ﬁnite diﬀerence schemes are often favored in cases where curvilinear multi-block grids can be generated, due to simpler coding and more eﬃcient use of computer resources. For aerodynamic applications where most of the surface of the aircraft is smooth, this methodology is especially suitable since i) curvilinear grids can be generated and ii) the resolution of large normal-to-surface gradients force the use of structured grids anyway. For very complicated geometries (such as close to landing gears), one can use hybrid methods (a combination of high-order ﬁnite diﬀerences and an unstructured method) as will be discussed below. Hybrid methods are also preferable in situations where waves propagate in free space after beeing generated by complicated geometrical features.
In this article, we will review the progress made towards towards stable highorder ﬁnite diﬀerence schemes for ﬂuid dynamics as well as other applications. To this end, we will brieﬂy explain the basic principles in a few examples. We will also discuss the SBP-SAT interpretation of other schemes and recent extensions of SBP-SAT schemes for time integration, non-linear theory and shock capturing.
The article is organized as follows. In Section 2, we present the theory for linear initial-boundary-value problem. We introduce the SBP-SAT concepts via a number of examples in Sections 3.1, 3.2 and 3.3. In Section 3.4 we discuss convergence rates and in Section 3.5 alternative ways to impose boundary conditions. In Section 3.6 and 3.7 we explain the SBP-SAT method in a 2-D example. In Section 3.8 we discuss aspects of the time evolution of the discrete system and in Section 3.9 we review results regarding dual consistency. In Section 4 we relate the SBP theory for ﬁnite diﬀerence schemes to other numerical methods. Section 5 contains a review of the various applications where SBP-SAT schemes have been used to obtain numerical approximations. Finally, we discuss some aspects of non-linear theory in Section 6.

2. Theory for Initial Boundary Value Problems
We begin by reviewing the general theory for Initial-Boundary-Value Problems (IBVP). Most of the material in this section can be found in [GKO95]. This sets the scene for the subsequent sections focusing on SBP-SAT schemes.

2.1. Preliminaries. Consider the initial-boundary-value problem

ut = P (x, t, ∂x)u + F, 0 ≤ x ≤ 1, t ≥ 0,

(1)

u(x, 0) = f (x),

L0(t, ∂x)u(0, t) = g0(t), L0(t, ∂x)u(1, t) = g1(t),

where u = (u1, .., um)T and P is a diﬀerential operator with smooth matrix coeﬃcients. L0 and L1 are diﬀerential operators deﬁning the boundary conditions. F = F (x, t) is a forcing function.

Deﬁnition 2.1. The IBVP (1) with F = g0 = g1 = 0 is well-posed, if for every f ∈ C∞ that vanishes in a neighborhood of x = 0, 1, it has a unique smooth solution

SUMMATION-BY-PARTS

3

that satisﬁes the estimate (2)

u(·, t) ≤ Keαct f

where K, αc are constants independent of f .
A problem is well-posed if it satisﬁes an estimate like (2). This require that appropriate boundary conditions are used which, along with the estimate, guarantees that a unique smooth solution exists. The extension to inhomogeneous boundary condition is possible via a transformation u˜ = u − Ψ where Ψ(x, 0) = f (x) and Ψ({0, 1}, t) = g0,1 such that u˜ satisﬁes (1) with homogeneous data (and a diﬀerent but smooth forcing function). However, to obtain Ψ, g0,1 is required to be differentiable in time. This requirement is not necessary if the problem is strongly well-posed as deﬁned below.

Deﬁnition 2.2. The IBVP (1) is strongly well-posed, if it is well-posed and

t

(3)

u(·, t) 2 ≤ K(t) f 2 + ( F (·, τ ) 2 + |g0(τ )|2 + |g1(τ )|)dτ

0

holds instead of (2). The function K(t) is bounded for every ﬁnite time and independent of F, g0, g1, f .

Remark We have tacitly assumed that boundary and initial data are compatible. Compatibility is necessary to ensure a smooth solution. E.g., for continuity one must require that g0(t) = f (0) and g1(t) = f (1). For higher continuity, derivatives of g0,1 and f must be related via the equation. (See [GKO95].)
Next, we turn to semi-discretizations of (1). Let xj = jh, j = 0, 1, ..., N where h = 1/N is the grid spacing. We deﬁne the grid functions fj = f (xj ) and Fj (t) = F (xj , t). With each grid point, we also associate a function (the approximate solution) vj(t). We will use the notion smooth grid function to denote a grid function being the projection of a smooth function. Furthermore, we use · h to denote a discrete L2-equivalent norm.
Then we approximate (1) by

(4)

(vj)t = Qj(xj , t)vj + Fj + Sj, j = 0, ..., N, t ≥ 0

vj(0) = fj

where Qj is the approximation of P at xj. Sj is the SAT term. Sj is zero except at a few points close to the boundary. The next deﬁnition is in analogy with Deﬁnition 2.1 above.

Deﬁnition 2.3. Consider (4) with F = g0 = g1 = 0. Let f be the projection of a C∞ function that vanish at the boundaries. The approximation is stable, if, for all h ≤ h0

(5)

v(t) h ≤ Keαdt f h

holds and K, αd are constants independent of f .

The same arguments as in the case of well-posedness can be used here to extend this notion to general inhomogeneous data in L2. To do so, g0 and g1 must be diﬀerentiable. The following deﬁnition rids this assumption.

Deﬁnition 2.4. The approximation (4) is strongly stable if it is stable and

(6)

v(t)

2 h

≤

K (t)

f

2 2

+

max
τ ∈[0,t]

F (τ )

2 h

+

max
τ ∈[0,t]

g0(τ )

2 h

+

max
τ ∈[0,t]

g1(τ )

2 h

holds instead of (5). K(t) is bounded for any ﬁnite t and independent of F, g0, g1, f .

4

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

Let u¯ be the projection of the exact solution u(x, t) on the grid, i.e., u¯i = u(xi, t). Then the (local) truncation error, T , is deﬁned by

(u¯j)t = Q(xj , t, D)u¯j + Fj + Sj + Tj, j = 0, ..., N, t ≥ 0

For an SBP-SAT scheme the truncation error takes the form,

(7)

T = (O(hr), ...O(hr), O(hp)...O(hp), O(hr)...O(hr))T

where r < p and the number of points with accuracy r is ﬁnite and conﬁned to the vicinity of the boundary. A scheme with r, p > 0 is termed consistent. The order of accuracy refers to the exponent in the truncation error. In (7) that would be r (if r < p) or in this case we may say (p, r) given that the structure of T is known. Moreover, the (solution) error is deﬁned as ej(t) = vj − u¯j and the convergence rate of the method is q if e 2 ≤ O(hq), where · 2 denotes the discrete L2 norm.
Although the deﬁnitions of (strong) well-posedness and (strong) stability are similar, the bounds in the corresponding estimates need not be the same. The following deﬁnition connects the growth rates of the continuous and semi-discrete solutions.

Deﬁnition 2.5. Assume that (1) is well-posed with αc in (2) and that the semidiscrete approximation is stable with αd in (5). If αd ≤ αc + O(h) for h ≤ h0 we say that the approximation is strictly stable.

In the above deﬁnitions, the schemes are semi-discrete, i.e., time is left continuous. Clearly, only fully discrete schemes are useful in practice. In [KW93], it was shown that semi-discrete stable schemes are, under certain conditions, stable when discretized in time using Runge-Kutta schemes. This problem was further studied in [LT98] with a particular focus on energy stable semi-discrete schemes (such as the scheme above). Recently it was shown in [NL13] how to extend the the SBP-SAT technique in space to the time-domain, where fully discrete sharp energy estimates can be obtained, see section 5.1.1 for more details.

3. Theory for SBP-SAT schemes In this section we present the SBP-SAT schemes by considering a few examples.

3.1. The advection equation. Consider

(8)

ut + aux = 0, 0 ≤ x ≤ 1, t > 0,

u(x, 0) = f (x),

u(0, t) = g0(t)

where f and g0 are initial and boundary data in L2 and a > 0. The semidiscretization of (8) is,

(9)

ut + aDu = P −1S, t > 0

u(0) = f

where u(t) = (u0(t), u1(t), ...uN (t))T and similarly f = (f (x0), ..., f (xN ))T . S = (S0, 0, ..., 0)T is the SAT enforcing the boundary condition (to be deﬁned later). D is a diﬀerence matrix.

Deﬁnition 3.1. D is a (p, r)-accurate ﬁrst-derivative SBP operator, if its trunca-

tion error is given by (7), D = P −1Q where Q + QT = B = diag(−1, 0, 0, ..., 0, 1),

P is symmetric positive deﬁnite and deﬁnes an L2-equivalent norm

v

2 P

= vT P v.

SUMMATION-BY-PARTS

5

The SBP operators for ﬁrst-derivative approximations associated with explicit central diﬀerence schemes are found in [Str94]. (Note that the entries of P scale as the stepsize h.) In [KS74] it was proven that if P is a diagonal matrix, the boundary closures can be at most of order r = τ when the interior stencils are of order p = 2τ . Such operators exist with order up to eight in the interior. If P is allowed to be non-diagonal with small blocks at the boundaries, the truncation error can be (τ, τ − 1). Furthermore, there exists SBP schemes for other than explicit ﬁnite diﬀerences. In [AC00, ACY00] and [CGA94], SBP-type boundary closures for compact interior discretizations were derived.
Let a+ = max(a, 0) = (|a| + a)/2. Stability of the advection equation was established in [CGA94] for SBP-SAT schemes.

Proposition 3.2. Let D be an SBP operator and S0 = σa+[P −1]0(u0 − g0). If σ < −1/2 the scheme (9) is strongly stable.

Proof.

u

2 t

=

uTt P u + uP ut

=

−auT (Q + QT )u + 2uS

gives

u

2 t

≤

au20 +

2aσu0(u0 − g0). If σ < −1/2 we obtain

u

2 t

≤

−

aσ2 (1+2σ)

g02

.

3.2. The advection-diﬀusion equation. Consider

(10)

ut + aux = (ǫux)x, 0 ≤ x ≤ 1, t > 0.

u(x, 0) = f (x)

au(0, t) − ǫux(0, t) = g0(t) ǫux(1, t) = g1(t)

where both a and ǫ are positive constants.

Remark Equation (10) is often used as a model problem for Navier-Stokes equations and with that interpretation, the left and right boundary conditions are prototypes of far-ﬁeld conditions based on the characteristic direction a.

To demonstrate that the problem is strongly well-posed, we apply the energy method.

u

2 t

+

2ǫ

ux

2 = (au2 − 2ǫuux)0 − (au2 − 2ǫuux)1

= a−1 (au − ǫux)2 − (ǫux)2 0 − a−1 (au − ǫux)2 − (ǫux)2 1 .

The construction above where the boundary terms are factored into clean squares was ﬁrst done in [Nor95] and later used in [NC99]. Keeping in mind that a is positive we obtain the ﬁnal strong estimate

(11)

u

2 t

+

2ǫ

ux

2 = a−1

g02 − (ǫux)20

− a−1

(au − ǫux)21 − g12

.

The most straightforward SBP-SAT discretization of this problem is

(12)

ut + aDu = D(ǫDu) + P −1S,

where

(S0)0 = σ0(au0 − ǫ(Du)0 − g0(t)), (S1)N = σ1(ǫ(Du)N − g1(t)).

Proposition 3.3. The scheme is strongly stable with σ0 = −1, σ1 = −1.

Proof. Use the energy method to obtain

u

2 t

+

2ǫ

Du

2 = a−1

g02 − (au − g0)20

− a−1

(au − g1)2N − g12

,

which is very similar to the continuous estimate (11) by considering the boundary conditions.

6

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

It should be remarked that the left boundary conditions is easily generalized to

βau(0, t) − ǫux(0, t) = g0(t),

β

>

1 2

.

The advection-diﬀusion equation in the SBP-SAT framework was ﬁrst studied in [CNG99], where the original proofs are found. Furthermore, the authors studied discretizations where the domain is subdivided in blocks, with possibly diﬀerent grid spacings and the stability conditions at the interfaces were derived.
We close this section by pointing out that other boundary conditions than those in (10) are possible. For instance in [SN08], Dirichlet boundary conditions were proven stable.

3.3. Interface treatment for the advection equation. To be able to generate a grid, it is often necessary to split the domain into several patches, each containing one smooth piece of the grid. (See the Section 1 and 5.1 for discussions on multiblock gridding.)
Here, we will demonstrate the interface procedure for the advection equation. Consider,

(13)

ut + aux = 0, −1 ≤ x ≤ 1, t > 0.

u(x, 0) = f (x)

au(−1, t) = ag−1(t)

where a > 0 is a constant. In Section 3.1, we showed that this problem is well-posed. (The change in domain size is not important for this conclusion.) Furthermore, a stable single-grid discretization, including the boundary condition, was proposed. Here, we split the domain at x = 0. We discretize the left piece, −1 ≤ x ≤ 0 with N + 1 points and the right, 0 ≤ x ≤ 1, with M + 1. We denote the grid spacings as hL = 1/M and hR = 1/N . The subscripts L and R will be used to signify operators in the two domains, e.g. DL = PL−1QL is the SBP diﬀerence operator in the left domain. We denote the discrete solution vectors as v in the left domain and u in the right. Note that both vN (t) and u0(t) represent approximations of u(0, t). Hence, the interface condition vN (t) = u0(t) will be enforced weakly in the scheme. We write the semi-discrete scheme as,

vt + aDLv = SL + S0L, t > 0

(14)

v(0) = fL,

ut + aDRu = S0R,

u(0) = fR.

Here SL is the SAT at x = −1 and corresponds to S in Proposition 3.2.

Proposition 3.4. Let DL, DR be SBP operators, [S0L]N = σL[PL−1]N (vN − u0)

and

[S0R]0

=

σR[PR−1]0(u0 − vN )

with

σL

≤

a 2

and

σR

= σL − a.

Then

the

scheme

(14) is strongly stable. Furthermore, the discretization is conservative in the sense

that it satisﬁes the weak form of the equation.

Proof. In the energy estimate we lump all the boundary terms at x = −1 in a term BTL and conclude that they are bounded due to Proposition 3.2. The energy estimate is obtained by multiplying the left scheme by vT PL and the right by uT PR and adding the results,
( v 2)t + ( u 2)t = −avN2 + au20 + 2vN σL(vN − u0) + 2u0σR(u0 − vN ) + BTL.
Using σL ≤ a/2 and σR = σL − a yields negative coupling terms. Hence, ( v 2)t + ( u 2)t ≤ BTL and we have a bound.

SUMMATION-BY-PARTS

7

Remark Note that the sign of a is immaterial to the interface part of this proof. The SAT terms at the interface are identical if a < 0

Next, we turn to conservation. Speciﬁcally, we show that the discrete scheme will, upon convergence, capture a weak solution. To this end, we introduce a smooth test function Φ that vanishes at the boundaries. The weak form of (13) is

(15)

1

t1

Φu dx|t0 −

(Φtu + Φxau) dxdt = 0

−1

0 −1

We denote its restriction to the grid as φL,R such that (φL,R)j(t) = Φ(xj, t) and (φL)N = (φR)0. We multiply the left part of the scheme (14) by φT PL, the right by φT PR and integrate in time to obtain,

(16)

φTL PLv|T0 + φTRPRu|T0 −
T
(φL)Tt PLu + (φR)Tt PRu + a(DLφL)T PLv + a(DRφR)T PRu dt = 0
0

where all the terms at the interface cancel out due to the conservation condition σR = σL − a. The relation (16) converges to (15) thanks to the L2 bound on the discrete solution.

The proof here is a special case of the one given in [CNG99]. (See also [GN11,
BN12a, EAN11].) In [CNG10] a general treatment of interfaces is presented. Non-
conforming interface procedures are derived in [NKG12, MC10]. For generalizations
to the multi-D Euler and Navier-Stokes equations, we refer to [NC99, NGvdWS09]. Non-linear analysis of interfaces can be found found in [FCN+13, SO¨ 13]

3.3.1. Second-derivative SBP operators. In (12) the second derivative is discretized by applying the ﬁrst derivative approximation twice. This leads to a wide stencil for the second derivative approximation. To increase the eﬀective accuracy, more compact stencils approximating the second derivative can be used. That is (assuming ǫ = constant), we would like a discretization

(17)

ut + aDu = ǫD2u + S

For this approximation to be SBP, certain conditions on D2 must be imposed. They were originally identiﬁed in [CNG99].

Deﬁnition 3.5. The second-derivative SBP matrix is deﬁned as D2 = P −1(−ST M + B)S with the following properties. Let u be a smooth grid function, then D2u = uxx + O(hr) and [Su]0,N = ux(x0,N ) + O(hr). Furthermore, M is symmetric positive deﬁnite with hnI ≤ M ≤ N hI where n, N are some constants.

Clearly, D2 = DD = P −1QP −1Q = P −1((P −1Q)T P − B)P −1Q is one such derivative with S = P −1Q and M = P . Others are found in [CNG99, MN04]. The D2 operators have similar accuracy constraints as the ﬁrst-derivative approximations. That is, with diagonal P the order is (2τ, τ ) and with non-diagonal P , the boundary accuracy can be raised to (τ, τ − 1).
To accommodate a varying ǫ while still keeping the stencil narrow, a new set of operators were derived in [Mat11]. Furthermore, in [MSS08] the accuracy and stability properties of various diﬀerent SBP second-derivative approximations were discussed.
The papers [AD97, AD99] are also SBP-SAT-like schemes for diﬀusive equations, but designed as an immersed boundary technique. That is, the domain is covered by a non-conforming Cartesian grid.

8

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

Remark An often expressed criticism against the use of DDu is that the stencils are wide and that they do not damp the highest frequency and therefore are prone to ”spurious oscillations”. That is only true for non-linear problem in the nonlinear regime for which the stability proofs do not hold. For linear problems such oscillations can not be generated by a stable scheme due to the bound in L2 ∩ L∞.
3.4. Convergence rates. The famous Lax-Richtmyer Theorem ([LR56]) states that a consistent scheme is convergent if it is stable. The theorem is very generally posed and does not require data (and consequently the solution) to be smooth. Therefore, no prediction of the convergence rate is possible. Our deﬁnitions of well-posedness requires solutions to be smooth, which in turn allows estimates of convergence rates. However, deriving convergence rates is a non-trivial task. Considering, (1) and its discretization (4), it is easy to see that the error e satisﬁes a scheme equivalent to (4) and the estimates (5) or (6) apply to e. Therefore, if r = p, the convergence rate is trivially p. If on the other hand r < p, and since the estimates are stated in L2, the convergence rate can be shown to be at least min(r + 1/2, p). This rate, however, is suboptimal.
In [Gus75] and [Gus81] the question of global convergence rate was studied. These papers are commonly cited to motivate a drop in the boundary accuracy according to the rule: ”For a (p, r)-accurate approximation with r < p, the convergence rate is r + 1”. This is true given that the scheme satisﬁes a determinant condition (or equivalently the Kreiss condition [GKO95]). However, it is not easy to prove that the Kreiss condition is satisﬁed for a realistic problem, and hence the analysis is often neglected. It would be desirable to infer the higher convergence rate at the boundary directly from the energy estimate. An eﬀort in this direction is found, in [ADG00] where parabolic problems were studied and using the energy method the rate min(r + 3/2, p) was proven. In practice, however, a convergence rate of min(r + 2, p) is observed.
In [SN06] the question of convergence rate was revisited with the intention to tie the convergence rate to energy estimates. For schemes approximating PDEs with principal part of order q and under certain accuracy conditions for the boundary conditions, it was shown that the convergence rate is min(r + q, p), if the numerical solution is bounded in L2 ∩ L∞. In a recent paper [Sv¨a13], L∞ bounds were derived using the energy estimates. Hence, any scheme satisfying an energy estimate is automatically bounded in L∞.
In (12), we approximate the second derivative by applying D twice. This makes the approximation of the second-derivative of order (r − 1, p), see [SN06]. Thanks to the bound in L2 ∩ L∞ ([Sv¨a13]) and [SN06], the global convergence rate is min(r + 1, p). The increased truncation error is exactly matched by the increased convergence rate for a parabolic problem. However, with the dedicated secondderivative approximations (17), the truncation error is (r, p) for both the hyperbolic and parabolic part and consequently the convergence rate is min(r + 2, p). For precise accuracy conditions on each operator, we refer to [SN06].
Remark Note that the results above hold even in the extreme case where the approximation of the PDE is inconsistent near the boundary. For instance, the biharmonic equation ut = −uxxxx was computed by applying the (4,3)-accurate ﬁrst-derivative 4 times. This rendered the approximation of uxxxx to be 0th-order accurate near the boundary. Still, a convergence rate of 4 was recorded in accordance with the theory.
3.5. Alternative boundary treatments. The SAT technique to enforce boundary conditions has been instrumental in the development of high-order ﬁnite diﬀerence schemes. It is relatively simple to implement and very versatile when proving

SUMMATION-BY-PARTS

9

stability. However, SATs are not the only way to enforce boundary conditions in conjunction with the SBP diﬀerence operators. The simplest, and most common, way to enforce boundary conditions for ﬁnite diﬀerence schemes is simply to overwrite the solution values on the boundary with boundary data every time step. We refer to this technique as the injection method. For simple PDEs, such as the advection equation, this may yield a stable approximation. However, already for linear systems in 1-D with waves travelling through the boundary in both directions this technique is doubtful. In [SM12] it was noted that injection was unstable for the Euler equations even at a supersonic inﬂow (where it would appear sensible to use since all characteristics are directed into the domain). For more information on the injection method, see [GKO95].
Another technique to enforce boundary conditions is the projection method. (See [Ols95a, Ols95b, GO96, GO98].) In this technique the solution is projected in each time step onto a subspace satisfying the boundary conditions and stability is proven with the energy method. For an extensive evaluation of the diﬀerent techniques to enforce boundary conditions, see [Mat03]. See also [Bod10] for a thorough evaluation of SAT boundary conditions in aeroacoustic applications; in [BN11] and also [SN08] further accuracy assessments of the SAT procedure are found.
Furthermore, the SAT method is a weak implementation of the boundary condition in the sense that the boundary points (u{0,N}) are unknowns and never explicitly set to the boundary value. Indeed, the boundary conditions are generally never satisﬁed exactly. The SAT acts as a forcing term driving this discrepancy to zero. On the other hand, both injection and projection are examples of strong enforcement of boundary conditions. Strong enforcement can sometimes be used for linear problems, if stability can be obtained. However, for the Euler equations, it is more subtle. It can be shown that to obtain a well-posed problem, the boundary conditions cannot be enforced strongly. (See [BF88]). This makes the SAT technique the preferred choice for non-linear problems.

3.6. Systems. We will introduce the use of SBP-SAT schemes in a multi-dimensional setting for a hyperbolic system of partial diﬀerential equations, which serves as a model for the 2-D Euler equations. For simplicity, we assume that the boundaries are all of far-ﬁeld type.

(18)

ut + Aux + Buy = 0, 0 < x, y < 1, 0 < t ≤ T

A+u(0, y, t) = A+gW (y, t)

A−u(1, y, t) = A−gE(y, t)

B+u(x, 0, t) = B+gS(x, t)

B−u(x, 1, t) = B−gN (x, t)

u(x, y, 0) = f (x, y)

Here, u is a column vector with m components and the matrices A, B are symmetric. Let A±, B± denote their positive and negative parts, i.e. XΛ+XT = A+ where Λ+ contains the positive eigenvalues and Λ+ + Λ− = Λ. We assume that gW,E,S,N , f are bounded in L2 ∩ L∞ and · denotes the L2-norm. The number of boundary
conditions are correct since only in-going characteristics are given at x, y = 0, 1.
First, we demonstrate well-posedness of (18) by deriving an energy estimate. That is we multiply (18) by uT and integrate by parts in space.

1 2

u

2 t

+

1 0

1
uT Auxdx dy +
0

1 0

1
uT Buydx dy = 0
0

10

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

or,

1

1

u

2 t

−

u(0, y, t)T Au(0, y, t)dy +

u(1, y, t)T Au(1, y, t)dy

0

0

1

1

− u(x, 0, t)T Bu(x, 0, t)dx + u(x, 1, t)T Bu(x, 1, t)dx = 0

0

0

In analogy with the SATs in the discrete case, we invoke the boundary conditions by adding

1

+2 u(0, y, t)T (A+u(0, y, t) − A+gW (y, t))dy (= 0)

0

1

(19)

−2 u(1, y, t)T (A−u(1, y, t) − A−gE(y, t))dy (= 0)

0

1
+2 u(x, 0, t)T (B+u(x, 0, t) − B+gS(x, t))dy (= 0)

0

1

−2 u(x, 1, t)T (B−u(x, 1, t) − B−gN (x, t))dy (= 0)

0

Remark We could have added a scaling parameter in front of every boundary term in (19). The derivation below would then give a range of permissible values of the parameters.

To reduce notation, we only write down the derivation for x = 0 terms.

1

1

u

2 t

−

u(0, y, t)T A−u(0, y, t)dy +

u(0, y, t)T A+u(0, y, t)dy

0

0

1

−2 u(0, y, t)T A+gW (y, t)dy + ... = 0

0

The terms associated with outﬂow at each boundary do not contribute to a growth and are dropped. (The A− term at x = 0 etc.) We obtain

(20)

1

1

u

2 t

≤

gW T A+gW dy − (u − gW )T A+(u − gW ) dy.

0

0

The ﬁrst term on the right-hand side is positive but bounded and represents the energy that is passed into the domain through the boundary. The second term vanishes thanks to the boundary condition. Integrating in time, from 0 to the ﬁnal time T , will give the desired bound on u(·, ·, T ) .

Remark No-penetration wall boundary conditions can be treated in a similar way and an energy estimates will follow, [SN08].

Remark If A and B are obtained from the linearized Euler equations, the above estimate demonstrates stability of smooth solutions to the non-linear Euler equations.

3.7. Spatial discretization. To deﬁne an SBP-SAT semi-discretization of (18),
we introduce the computational grid, xi = ihx, i ∈ {0, 1, 2, ..., N } and yj = jhy, j ∈ {0, 1, 2, ..., M } where hx, hy > 0 are the grid spacings. We will also need the vectors e0 = (1, 0, 0..., 0)T and eN = (0, ..., 0, 1)T , and the matrices E0 = e0eT0 (1 in the upper-right corner and 0 elsewhere) and EN = eN eTN (1 in the lower-left corner and 0 elsewhere). With this notation we have Q + QT = −E0 + EN .

SUMMATION-BY-PARTS

11

Next we extend the SBP operators to two dimensions by the use of Kronecker products. The Kronecker product of two matrices A, B is deﬁned as

 a11B . . . a1N B 

(21)

A⊗B = 

...

...

. 

aN1B . . . aNN B

It satisﬁes A ⊗ B + C ⊗ D = (A + C) ⊗ (B + D) and (A ⊗ B)(C ⊗ D) = (AC ⊗ BD)
assuming that both A and C have the same dimensions and so has B and D. Furthermore, (A ⊗ B)T = (AT ⊗ BT ) and, if A, B are invertible then (A ⊗ B)−1 = (A−1 ⊗ B−1).
Let IM denote an M -by-M identity matrix and let vkij denote the numerical approximation of uk(xi, yj), i.e., the kth component of u at xi, yj. We introduce:

v = (v111, ...vm11, v121...vm21, ...vmN1, v112, ...vm12, ...)T

P = (Py ⊗ Px ⊗ Im)

P deﬁnes a 2D l2-norm by v 2 = vT Pv. Note that (IM ⊗ DN ⊗ Im)v calculates the x-derivative approximation in the entire domain. Also, (IM ⊗ E0 ⊗ Im)v extracts the boundary points at x = 0, i.e., it sets all but the boundary points to 0 in v. With these observations, we write the scheme as

(22)

vt + (IM ⊗ Dx ⊗ A)v + (DM ⊗ Ix ⊗ B)v = SAT

where (23)
(24)

SAT = − (IM ⊗ Px−1E0 ⊗ A+)(v − g) + (IM ⊗ Px−1EN ⊗ A−)(v − g) − (Py−1E0 ⊗ IN ⊗ B+)(v − g) + (Py−1EM ⊗ IN ⊗ B−)(v − g)

g has the same structure as v with gW,E,S,N injected on the appropriate boundary positions and it is 0 elsewhere.
To reduce notation below, we single out one boundary

(25)

SAT = −(IM ⊗ Px−1E0 ⊗ A+)(v − g) + SAT

Next, we derive a bound on v. That is, we multiply (22) from the left by vT P.

vT Pvt + vT (Py ⊗ Qx ⊗ A)v + vT (Qy ⊗ Px ⊗ B)v = −vT (Py ⊗ E0 ⊗ A+)(v − g) + vT PSAT

(Notice that the ﬁrst term on the right-hand side corresponds exactly to the ﬁrst in (19).) Adding the transposed equation, and using QN + QTN = −E0 + EN , yields,
(vT Pv)t + vT (Py ⊗ −E0 + EN ⊗ A)v + vT (−E0 + EM ⊗ IN ⊗ B)v =
−2vT (Py ⊗ E0 ⊗ A+)(v − g) + 2vT PSAT

To reduce notation, we set all boundary terms except those at x = 0 to 0. (The others are dealt with in the same way.) We get,

v

2 t

= vT (Py ⊗ E0 ⊗ (A+ + A−)v − 2vT (Py ⊗ E0 ⊗ A+)(v − g)

v

2 t

≤ −vT (Py ⊗ E0 ⊗ A+)v + 2v(Py ⊗ E0 ⊗ A+)g

Denote M = (Py ⊗ E0 ⊗ A+) and note that the matrix is bounded and positive semi-deﬁnite. We obtain

(26)

v

2 t

≤ −vT M v + 2vM g = gT M g − (v − g)T M (v − g),

12

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

which results in a bounded growth and hence the scheme is stable. The right-hand side (26) is a direct analog of (20). The diﬀerence here is that v is not necessarily equal to g and hence a small dissipation is added at the boundary.
The SAT (25) could be stated more generally as,

(27)

SAT = −(Iy ⊗ Px−1E0 ⊗ A˜)(v − g) + SAT

where we have shown stability for A˜ = A+ which represents a far-ﬁeld boundary. By choosing A˜ diﬀerently, we may enforce a no-penetration condition and prove stability in a similar way. Interfaces between diﬀerent grids also ﬁt in this framework; A˜ is chosen to obtain an energy estimate and g would contain data from the neighboring grid block. We stress that an interface is not a boundary in its ordinary meaning. Therefore, there is no constraint with regards to the number of conditions given (although too few, meaning less than in-going characteristics, will lead to an unstable approximation).
For more information on boundary procedures and generalizations to the NavierStokes equations, see [NGvdWS09, SN08, SCN07, GN11].

3.8. Strict stability. So far we have discussed standard energy estimates which produce L2-bounds on the solution at any ﬁnal time T . Such estimates ensure stability in the sense that on any compact domain (in space and time) the solution converges as the grid is reﬁned. For a particular grid resolution, however, the time evolution of the energy may diﬀer between the numerical and true solution. In Deﬁnition 2.5, strict stability was deﬁned to ensure a numerically correct time evolution of the discrete solution; the time evolution of the norm should converge with h as the grid is reﬁned. Sometimes it is desirable to take this requirement one step further by demanding that αd = αc − |ch| such that the discrete scheme is more dissipative than the PDE itself. For constant coeﬃcient PDEs on a Cartesian domain, this is generally not an issue since strictly stable approximations are usually obtained. For variable coeﬃcient PDEs, this is not the case.

3.8.1. Splitting. Consider a variable coeﬃcient scalar PDE

(28)

ut + (a(x)u)x = 0, 0 < x < 1

augmented with appropriate boundary and initial conditions. We assume that a(x) is smooth such that maxx |ax| is bounded. For smooth solutions we can rewrite (28) on skew-symmetric form

(29)

ut

+

1 2

(au)x

+

1 2

aux

+

1 2

ax

u

=

0.

The equations (28) and (29) are of course equivalent. We assume a(x) > 0 and give

a boundary condition to the left u(xL, t) = 0. Applying the energy method results in

xR

u

2 t

=

−a(xR)u2R

−

axu2 dx.

xL

The last term need not be negative since we have not assumed anything about the sign of ax. Hence, the energy may grow. However, it is bounded since,

(30)

u

2 t

≤

mxax(|ax|)

u

2.

We now have two choices of discretizations. The ﬁrst is

ut + P −1Q(Au) = SAT.

SUMMATION-BY-PARTS

13

We omit further details and conclude that it is possible to prove stability by obtaining an energy estimate. (See [MS10] for estimates of variable coeﬃcient problems.) If, on the other hand, we choose the scheme

ut

+

1 2

D(Au)

+

1 2

ADu

+

1 2

Axu

=

0

where the entries of the matrix Ax are computed using P −1Qa ([a]j = a(xj )), we obtain the more accurate growth rate corresponding to (30). For more details on splitted schemes and strict stability, see [Nor06, Nor07]. Computational results demonstrating the eﬀect of splitting and diagonal norms (P diagonal) are found in [KDN13, KDN12] where the linear elastic wave equation is combined with nonlinear interface conditions (friction laws for faults) in earthquake modeling.
Similar splitting techniques can be exploited to obtain energy-style estimates for non-linear problems. Such splittings are generalizations of the well-known ”1/3”trick for Burgers’ equation, ut + (u2/2)x = 0. For smooth solutions this can be written as

ut

+

1 3

u2x

+

1 3

uux

=

0,

for which an estimate of u 2 is readily obtained. The key is homogeneity of the ﬂux as observed in [OO94], where estimates for the Euler equations were derived.
This was further developed in [GO96] and [SM12]. A ﬂux-splitting technique was also analyzed as a mean to stabilize non-linear equations in [FCN+13].

3.8.2. Coordinate transformations. Next, we will discuss the eﬀect of coordinate transformations on the time evolution of the semi-discrete problem. Consider

(31)

ut + ux = 0, u(0, t) = 0, 0 ≤ x ≤ 1

The energy method gives the following estimate: u 2 = −u2(1, t). The decay of the energy is precisely the energy disappearing through the right outﬂow boundary. If we discretize this problem with a constant grid size h, we get

ut

+

P

−1Qu

=

−

1 2

[P

−1]0(u0

−

0)

which results in the estimate

u

2 t

= −u2N ,

i.e.,

exactly

the

same

as

the

continuous

estimate. The discretization is strictly stable.

Next introduce a (non-singular) coordinate transformation, x = x(ξ) (e.g. a

stretched grid). Then, we obtain the equivalent equations

(32)

ut + uξξx = 0, u(0, t) = 0,

or,

(33)

(ξx)−1ut + uξ = 0, u(0, t) = 0

The problem (32) is now a variable coeﬃcient problem. We may discretize it as

(34)

ut

+

AP

−1Qu

=

−

1 2

[A11P

−1]0(u0

−

0)

where Aii = ξx(xi). To prove stability, we may freeze the coeﬃcients and the proof is essentially the same as for ut + Du = SAT above. However, when freezing the coeﬃcients we neglect terms that have an eﬀect, albeit bounded, on the energy

estimate. Hence, this discretization may not be strictly stable.

Another option is to consider (33) and note that due to the non-singularity of

x(ξ), A is positive deﬁnite. From this we may introduce a weighted L2 norm,

v

2 ξ

= vT A−1P v

for P

diagonal.

In this norm, we obtain the

growth (

v

2ξ )t

=

−u2N for (34), which is the same as for (31) and hence the discretization is strictly

14

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

stable. It was proven in ([Sv¨a04]) that the condition that P is diagonal can not be circumvented.
This implies that while any SBP operator used to compute (34) is stable, it is only those with P diagonal that can be proven to have the same time evolution as the continuous problem. For instance, if a block-norm SBP operator is used there may be small positive eigenvalues that cause a growth. This has been observed for the Euler equations on curvilinear meshes, where the block-norm operators require somewhat more artiﬁcial diﬀusion to converge to steady-state than their diagonalnorm counterparts, [SMN05]. A similar phenomena was observed in [KDN13].

3.9. Dual consistency and superconvergence of functionals. In many cases, accurate solutions to the equations themselves might not be the primary target for a calculation. Typically, functionals computed from the solution, such as the lift and drag coeﬃcients on a body in a ﬂuid, potential or kinetic energy, or probabilities are of equal or even larger interest. The importance of duality in the context of functionals was realized in adaptive mesh reﬁnement, error analysis and optimal design problems, which has made the study of duality somewhat restricted to unstructured methods such as ﬁnite element (FEM), discontinuous Galerkin (DG) and ﬁnite volume (FVM) methods.
Recently, however, it was shown in [HZ11, Hic12, HZ13] that the adjoint equations can be used for ﬁnite diﬀerence (FD) methods to raise the order of accuracy of linear functionals of numerical solutions generated by SBP-SAT schemes. In general, the truncation error of SBP-SAT discretization with diagonal norms is of order 2p in the interior and p at the boundary, which results in a solution accuracy of p + 1. (See [SN06].) It was shown in [HZ11], and later extended in [BN12b],[HZ14], that linear integral functionals from a diagonal norm dual-consistent SBP-SAT discretization retains the full accuracy of 2p. That is, the accuracy of the functional is higher than the solution itself. This superconvergent behavior was previously observed for FEM and DG methods but it had not been previously proven for ﬁnite diﬀerence schemes.
It is important to note that dual consistency is a matter of choosing the coefﬁcients in the SATs and it does not increase the computational complexity. Superconvergence of linear integral functionals hence comes for free. Free superconvergence is an attractive property of a dual consistent SBP-SAT discretization but the duality concept can also be used to construct new boundary conditions for the continuous problem. Research in this direction was initiated in [BN13b], and the procedure is under development for both the Euler and Navier-Stokes equations [BN13a].

3.10. SBP-SAT in time. In [NL13] the SBP-SAT technique in space was extended to the time-domain. To present the main features of the methodology, we consider the simplest possible ﬁrst order initial-value problem

(35)

ut = λu,

with initial condition u(0) = f and 0 ≤ t ≤ T . Let λ be a scalar complex constant representing an energy stable spatial discretization of an IBVP. Hence, we assume that Re(λ) < 0.
Remark Note that the eigenvalues of diﬀerence approximations may be complex motivating our choice of λ. Furthermore, the size of the maximal eigenvalue of a diﬀerence approximation of a hyperbolic problem grows as 1/h and for a parabolic problem as 1/h2.

SUMMATION-BY-PARTS

15

The energy method (multiplying with the complex conjugated solution and integrating over the domain) applied to (35) yields

(36)

|u(T )|2 − 2Re(λ)||u||2 = |f |2,

where ||u||2 =

T 0

|u|2dt.

Since

Re(λ)

< 0,

the

solution

at

the

ﬁnal

time

is

bounded

in terms of the initial data and we also obtain a bound on the (temporal) norm of

the solution.

An SBP-SAT approximation of (35) reads

(37)

P −1QU = λU + P −1(σ(U0 − f ))e0.

The vector U contains the numerical approximation of u at all grid points in time. The matrices P, Q form the diﬀerentiation matrix with the standard SBP properties given in Deﬁnition 3.1 The penalty term on the right-hand-side of (37) enforces the initial condition weakly using the SAT technique at grid point zero using the unit vector e0 = (1, 0, ..., 0, 0)T .
Remark The penalty term in (37) forces the discrete solution towards the initial data, i.e. U0 = f in general, but it is close.
The discrete energy method applied to (37) (multiplying from the left with U ∗P , using the SBP properties and making the choice σ = −1 leads to

(38)

|UN |2 − 2Re(λ)||U ||2P = |f |2 − |U0 − f |2.

The choice σ = −1 also makes the discretization dual consistent [HZ11],[BN12b],

[BN13b],[HZ14]. By comparing the continuous estimate (36) with (38) we see that

the discrete bound is slightly more strict than the continuous counterpart due to the term −|U0 − f |2 (which goes to zero with increasing accuracy).

Remark Note that the estimate (38) is independent of the size of the time-step,

i.e. the method is unconditionally stable.

Remark Sharp estimates like (38) can hardly be obtained using conventional local methods where only a few time levels are involved. One can argue, although no proof exist, that it can be done only with global methods.

In [NL13] it is shown that the new SBP-SAT method in time is high order accurate, unconditionally stable and together with energy stable semi-discrete approximations, it generates optimal fully discrete energy estimates. In particular, for energy stable multi-dimensional system problems such as the Maxwells equations, the elastic wave equations and the linearized Euler and Navier-Stokes, fully discrete energy estimates for high order approximations can be obtained in an almost automatic way.
In [LN13], it was shown how the SBP-SAT technique for time integration, originally formulated as a global method, can be used with ﬂexibility as a one-step multi-stage method with a variable number of stages proportional to the order of the scheme, without loss of accuracy compared to the global formulation. This fact makes the SBP-SAT method highly competive, easier to program and signiﬁcantly reduce the storage requirements. Classical stability results, including A- stability, L-stability and B-stability could also be proven using the energy method. In addition is was shown that non-linear stability holds for diagonal norm operators when applied to energy stable initial value problems.

4. SBP and other numerical schemes
In this section, we will digress from the high-order ﬁnite diﬀerence scheme discussed so far. As stated in Deﬁnition 3.1, the SBP property is not tied to diﬀerence schemes but to the properties of certain matrices. A numerical derivative resulting from any numerical discretization technique can be expressed on matrix form and

16

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

hence, it may or may not have the SBP property. If it has, we can employ the SAT framework to enforce boundary conditions and eﬀectively reuse all the stability theory stemming from the research on high-order SBP ﬁnite diﬀerence schemes.

4.1. Finite volume schemes. In [NFAE03] it was shown that the unstructured node-centered ﬁnite volume ﬁrst-derivative approximation is on SBP form. In [NB01] also the structured cell-centered ﬁnite volume method is modiﬁed and extended to SBP form. In [SN04] it was shown that a common Laplacian approximation is also on SBP form and SBP consistent artiﬁcial diﬀusion was proposed in [SGN06]. The entire framework for imposing boundary conditions with SAT terms was successfully used in [SSHM10] for ﬁnite volume schemes. The weak enforcement of boundary conditions created the possibility of hybrid ﬂow solvers where ﬁnite volume and ﬁnite diﬀerence schemes are used in diﬀerent domains. This greatly simpliﬁes grid generation since complicated geometries can be wrapped in unstructured grids and in the bulk of the domain structured grids and the more accurate and eﬃcient ﬁnite diﬀerence schemes are used. Hybrid schemes schemes have been developed in [Gon06, GN07, NHS+09].

4.2. The discontinuous Galerkin method. The SBP-SAT method can also be

related the discontinuous Galerkin (dG) method. In fact, one may argue that it is

the same methodology, only in a diﬀerent technical setting. Consider the advection equation (8). We expand the solution in a polynomial u = LT (x) α(t) where L(x) = (φ0, φ1, ... φN )T and α(t) = (α0, α1, ... αN )T . By inserting u into (8),
multiplying each row with the basis functions φi and integrating over the domain
we ﬁnd

1

1

(39)

LLT dxαt + a LLTx dxα = 0, ⇒ P αt + aQα = 0.

0

0

The matrices P and Q satisfy the SBP requirements in Deﬁnition 3.1 if one chooses

φi to be the ith Lagrange polynomial. For more details, see [Nor06].
How about imposing boundary conditions ? To see the similarities we replace Q in (39) by −QT + B (integration-by-parts with SBP operators). We obtain the modiﬁed equation P αt + aBα − aQT α = 0. Next we use a common dG procedure
and replace”what we have” (α0) with ”what we want” (g(t)) and integrate back by replacing QT with −Q + B. The ﬁnal result is

(40)

P αt + aQα = −a(α0 − g(t))e0

where eT0 = (1, 0...0). Clearly (40) includes a weak SAT term for the boundary condition and consequently the dG scheme presented above is on SBP-SAT form.

For more details on relations to dG schemes, see [Gas13]. Other striking similarities

with SBP-SAT schemes can be found in the so called ﬂux reconstruction schemes,

where the penalty terms are constructed using specially designed polynomials. For

a description of that method and similar ones, see [WCVJ13].

5. Applications
5.1. The multi-dimensional Navier-Stokes equations. In this section, we will brieﬂy sketch the procedure for the Navier-Stokes equations and we refer to the references for more details. In general, the domain is not Cartesian. To be able to use ﬁnite diﬀerence schemes, it must be possible to subdivide the computational domain in several non-overlapping patches that cover the entire domain. In each subpatch, a curvilinear grid is constructed. The grid must be a suﬃciently smooth to allow for a smooth transformation to a Cartesian mesh. That is, we assume that there are coordinates 0 ≤ ξ, η ≤ 1 and a transformation x = x(ξ, η), y = y(ξ, η) to the physical space in each patch. Grid lines should be continuous across an

SUMMATION-BY-PARTS

17

interface but their derivatives need not. We remark that grid lines are allowed to end at an interface, i.e., the resolution of two neighboring blocks may be diﬀerent, see [MC10] and also [NKG12].
In each patch, the transformed system takes the form,

(41)

J ut + (A˜1u)ξ + (A˜2u)η = ǫ(F˜ξvξ + F˜ηvη), 0 < ξ, η < 1, t > 0

F˜vξ = B˜11uξ + B˜12uη

F˜vη = B˜21uξ + B˜22uη

augmented with suitable boundary, interface and initial conditions. The process to obtain an energy estimate begins with a linearization followed
by a symmetrization, see [AG81]. To state the boundary conditions, the matrices, A˜±1,2, are needed. The diagonalization of these matrices were carried out in [PC81]. Derivation of energy estimates for the Euler and Navier-Stokes equations with various types of boundary and interface are found in [NC99, NC01] where the speciﬁc diﬃculties associated with these equations are explained. Furthermore, in [SCN07] far-ﬁeld boundary conditions are derived and particular attention given to the diﬃcult case of subsonic outﬂow. In [SN08] solid wall and in [BN11] Robin wall boundary conditions are derived. Also, recall that for curvilinear grids, SBP operators with diagonal P , have favorable stability properties as discussed in Section 3.8.2 and [Sv¨a04]. The speciﬁc treatment of grid block interfaces is addressed in [NGvdWS09].
The diﬀusive terms in the Navier-Stokes equations may be discretized in a few diﬀerent ways. In (41) above, they are treated as ﬂux terms and hence they are approximated by repeated applications of the ﬁrst-derivative approximation. As discussed in section 3.3.1 it may be advantageous from an accuracy point of view to use dedicated second-derivative approximations that use narrow stencils. To this end, we rewrite (41) as

J ut + (A˜1u)ξ + (A˜2u)η = ǫ(B˜11uξξ + B˜12uηξ + B˜21uξη + B˜22uηη).

This derivation assumes that the matrices B˜ij are constant. We can then use the second-derivative operator given by D2 = P −1(−ST M S + BS) on uξξ and uηη. (Recall that BS is a matrix with the ﬁrst and last rows non-zero and approximating
the ﬁrst-derivative.) The cross-derivative terms are still approximated with the ﬁrst-derivative, P −1Q, and we note that such an approximation is not ”wide” since
they operate in diﬀerent directions. To obtain an energy estimate, the D2 and D1 operators must relate in a certain way. Roughly speaking, the D2 operator must be more diﬀusive than D1D1. Such D2 operators are termed ”compatible” and can be found in [MSS08]. In the case when the diﬀusion coeﬃcient is not constant,
D2 must be a function of the diﬀusion coeﬃcient. (C.f. (ǫux)x ≈ (ǫi+1/2(ui+1 − ui) − ǫi−1/2(ui − ui−1))/h2 in the interior of the scalar case with second-order approximation.) SBP versions of such approximations are found in [Mat11].

Remark It should be mentioned here that the use of compact second derivatives complicates the programming since one have to keep track on whether one deals with a clean second derivative or a mixed types of derivative. Next one have to build the ﬂux with these diﬀerent components. With the use of only ﬁrst derivative operators, no such bookkeeping is necessary, one simply build the complete ﬂux, and diﬀerentiate.

Many times it is necessary to damp small oscillations. This is usually done by adding artiﬁcial dissipation to the scheme. For high-order schemes such diﬀusion terms are even order diﬀerences scaled such that they do not to reduce the order of accuracy of the scheme and maintain a grid independent damping on the highest

18

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

modes. For instance, a fourth-order undivided second derivative is appropriate to add to a 4th-order accurate scheme. However, the artiﬁcial diﬀusion terms must be augmented with boundary closures that preserve the energy estimates. Without appropriate boundary closures, the artiﬁcial dissipation terms may destabilize the scheme. SBP conforming artiﬁcial diﬀusion operators can be found in [MSN04].

5.1.1. Other time integration techniques and computational aspects. To develop a code using the SBP-SAT schemes necessitates a number of choices apart from the spatial discretization scheme. For time-dependent problems, a time-discretization scheme must be chosen. In many of the computational results obtained in the references listed here a low-storage 3rd-order Runge-Kutta method has been used, [KCL00]. However, the time-step constraint for an explicit method is often severe. If the stiﬀness is artiﬁcial, i.e., not reﬂecting time-scales in the actual solution, an implicit scheme may be more eﬀective. We mention implicit-explicit Additive Runge-Kutta methods, [KC03], as a viable option where the cost of implicit schemes may be kept to a minimum by using implicit schemes in stiﬀ regions only. An eﬀort in this direction is found in [SSHM10] in the context of SBP ﬁnite-volume schemes. (See also [BCVK02, CKB+05].)
As a step towards eﬀective implicit schemes, eﬃcient steady-state solvers have been developed. While there are well-known and reliable libraries available for both linear and non-linear solvers there are still obstacles to overcome to use them in practice. One is to obtain the Jacobian matrix of the scheme. Already for second-order schemes it is cumbersome to derive and program exact Jacobians. For high-order methods even more so. Another option is approximate the Jacobian numerically. However, as the accuracy is increased by reducing the step size in the diﬀerence procedure, cancellation errors increase. To ﬁnd the optimal value is not easy and the optimal value may not be accurate enough for fast convergence to steady state. In [vdWS12], dual numbers are used (see [FJAvW11]) to approximate the Jacobian numerically. This is an algebraic trick that rids the problem with cancellation and fast convergence is obtained.
Finally, we return to the spatial discretizations, which naturally has an eﬀect on the convergence rate to steady state. It has already been mentioned that the use of diagonal-norm schemes (P diagonal) increases the convergence speed and allows for less diﬀusion in the simulation, [SMN05]. This has also been noted in [KDN13, KDN12], where it is further observed that splitting the convective terms has a favorable impact on convergence speed. Furthermore, the weak imposition of wall boundary conditions, i.e., the SAT procedure, increases the speed to steady state when compared to injecting the boundary conditions at the wall, [NEE12].

5.1.2. Computational Eﬃciency. We have to a great extent discussed stability proofs of SBP-SAT schemes. This is the foundation to ensure robust and highly accurate simulations. Next, we will review the literature regarding the eﬃciency of SBP schemes. In [SMN05, MSCN07] steady and unsteady airfoil computations governed by the Euler equations are found. In particular it is demonstrated that low order schemes fail (on reasonable grids) to propagate ﬂow structures over long distances. Airfoil computations (NACA 0012) with Navier-Stokes equations are found in [SLN10] and a subtle case with M a = 0.5 and Re = 5000 was successfully computed. The solution displays back ﬂow but the wake remain stable. Inaccurate simulations would produce an unstable wake. More simulations with the Navier-Stokes equations are found in [SCN07, SN08] for a ﬂow around a cylinder. The focus in these articles is on the robustness and accuracy of far-ﬁeld and wall boundary conditions. It is demonstrated that accurate shedding frequencies and separation points are obtained on course grids when the order of accuracy is suﬃciently high. Two

SUMMATION-BY-PARTS

19

other evaluations of the SAT technique to enforce boundary conditions are found in [Bod09, Bod10]. Here, the focus lies on aeroacoustic applications where minimal reﬂections from the boundaries are of the essence. As previously mentioned, the high-order SBP-SAT schemes was shown to be very competitive in terms of the use of computer resources, compared with other high-order techniques, see [vdWS12].
In [SSHM10], an unstructured ﬁnite volume SBP scheme was implemented for the Euler and Navier-Stokes equations. Implicit-explicit time integration were employed and the implicit part adaptively associated with stiﬀ regions in the problem. This procedure greatly reduced the computational cost and the eﬃciency was demonstrated on a number of test cases including LES around a cylinder.

5.2. Other applications. So far, we have mainly discussed SBP-SAT schemes applied to the compressible Navier-Stokes equations. Here, we will brieﬂy summarize other areas where SBP-SAT schemes have been successfully used. In [SM11], SBP operators were used to solve the compressible Euler equations with a stiﬀ source term modeling a reaction equation. Special boundary conditions that automatically produced divergence free solutions were derived and implemented in stable SBP schemes for the incompressible Navier-Stokes equations in [NMS07]. The incompressible Navier-Stokes equations were also considered in [HMIM07] where an accuracy evaluation of SBP-SAT ﬁnite volume schemes was the main objective.
Some applications in aeroacoustics, using the linearized Euler equations, are found in [Mu¨l08, MY02] and with focus on human phonation in [MY09]. Fluidstructure interaction problems were studied in [NE10], and it was shown that the weak coupling in SBP-SAT schemes was the coupling procedure that lead to the best result for this very sensitive problem. In [LN10] conjugate heat transfer between a solid and a ﬂuid was analyzed using SBP schemes. Moreover, in [NB13] the conjugate heat transfer was analyzed further by using the Navier-Stokes equations only (no heat equation was involved) but in combination with a special interface conditions. The resulting well-posed problem was implemented using an SBP-SAT scheme.
Another area particularly suited for high-order SBP schemes is wave propagation. In [MN06, MHI08, MHI09] the wave equation on second-order form was analyzed; boundary conditions as well as interface conditions between media with diﬀerent propagation speeds were derived and SBP-SAT schemes for solving the problem proposed.
In numerical relativity, the Einstein equations constitute a set of equations similar to the second-order wave equation and much eﬀort has gone into designing effective SBP-SAT schemes. A sample of articles in this ﬁeld is: [SDDT06, DDS+07, DBD+06, ZST08, MP10]. Furthermore, the SBP techniques have been applied to Maxwell’s equations of electromagnetism ([NG03]) and the magnetic induction equations ([KMRS09, KMRS12, MS10]). Recently, they have also been used in wave propagation problems in geophysics [KDN13, KDN12] and medicine [AN13].
The SBP-SAT methodology has also been used in the emerging area of Uncertainty Quantiﬁcation (UQ), which combines computational mathematics and mathematical statistics. Unlike the conventional approach using Monte Carlo simulations that require a large number of runs to obtain reliable statistics (and hence calls for extremely fast solvers), UQ is designed to solve directly for the solution and its associated statistics. Hence, UQ is potentially very eﬀective and SBP-SAT schemes produce accurate numerical solutions to the resulting PDE system. Problems in UQ solved the SBP-SAT technique are found in [PIN09, PNI10, PDN13, PIN13, PIN14].

20

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

6. Non-linear stability
The SBP-SAT schemes are often used to approximate non-linear PDEs, such as the Euler and Navier-Stokes equations. Before discussing non-linear stability, we will make a brief justiﬁcation for the linear stability proofs obtained with SBPSAT schemes. The Navier-Stokes equations can be linearized around a smooth solution (assuming it exists) which gives a variable coeﬃcient linear problem in the perturbed variables. If the solution of this problem is proven to be bounded then the non-linear problem is stable with respect to small perturbations. A similar argument is then made for a numerical discretization of the non-linear equations, [Str64]. The simpliﬁcation process is often taken one step further and the ”constant coeﬃcient” problem is analyzed. This means that the variable coeﬃcients in the PDE are ”frozen” prior to the analysis. This approach is justiﬁable, at least if the energy method is used to prove stability, see [KL89, MS10]. (If other techniques are used to prove stability the situation is much more subtle, see [Mic83, Mic87, Wad90, SW88]). This line of argument justiﬁes the convergence observed in subsonic regimes where non-linear eﬀects are negligible. Mild non-linear eﬀects can be kept at bay by the use of high-order diﬀusion terms, [MSN04, DDS+07].
Strong non-linear eﬀects in ﬂuid mechanics are often associated with shocks and the generation of non-smooth solutions for which the above linear arguments do not hold. However, it is possible to derive L2 bounds without prior linearization. Unfortunately, an L2 bound is not enough to imply convergence for a non-linear PDE. Nevertheless, stability of some sort is required and many diﬀerent paths have been taken.
A popular choice of ﬁnite diﬀerence schemes are WENO type schemes, originally derived in [JS96], that produce approximate solutions with sharp shock resolution and formally high order of accuracy. However, even linear stability without boundary conditions for such schemes are a subtle matter, [WS07, MMR11]. In two recent papers, the concepts of WENO and SBP were merged yielding linearly stable schemes for initial-boundary-value problems, [YC09, FCYF11].
Another approach was taken in [GO96, GO98] were they used a splitting of the ﬂux function based on entropy arguments (see [OO94]). Under certain assumptions, including a bounded wave speed, non-linear energy estimates were obtained. The drawback of this approach is that it usually introduces a non-conservative term which may cause an incorrect shock location. However, for scalar problems it was shown in [FCN+13] that most splitted schemes are in fact conservative, and hence the non-conservative terms will not cause a drift in the shock location.
A successful way to prove non-linear stability for conservation laws is reviewed in [Tad03]. The idea is to build a scheme that prevents the entropy from growing unboundedly and consequently they are termed entropy stable. The theory was developed for Cauchy problems and a global entropy estimate obtained using summation-by-parts (without boundaries). In [SM12], second-order SBP schemes were put in the entropy stability framework and non-linear stability proven for initial-boundary-value conservation laws, including the Euler equations. The theory was extended to 3rd-order accurate SBP schemes in [Sv¨a12]. Non-linearly entropy stable boundary conditions were revisited in [SO¨ 13] and here stability proofs for far-ﬁeld, solid wall and interfaces were given for the Euler equations. In the WENO framework, entropy stable schemes have been derived in [FC13] with SBP boundary closures such that the boundary conditions in [SO¨ 13, SM12] are readily applicable.

SUMMATION-BY-PARTS

21

7. Summary and Outlook
SBP-SAT schemes have reached a mature state where it is relatively straightforward to apply the method to new problems, as demonstrated by the vastly diﬀerent applications addressed to date. The mathematical foundation of SBP-SAT schemes, with the most important property being that convergence can be proved for linear or linearized problems, gives credibility to the numerical simulations that ad hoc methods can not oﬀer. The possibility to evaluate the size of numerical errors should be of paramount importance in both engineering and scientiﬁc simulations.
The key to the stability proofs is of course the boundary treatments of the SBP-SAT schemes. This is the feature that also makes the methodology versatile. It allows for stable and accurate couplings with other type of methods (hybrid schemes) and it is also the key to connect diﬀerent models in a stable manner. (E.g. ﬂuid-structure interaction.)
The remarks above outlines three paths for future research:
• Applying SBP-SAT to new problems which to a large extent involves analysis of boundary conditions.
• Connecting diﬀerent models which calls for analysis of interface conditions and the subsequent analysis to derive stable SAT connections
• Derivation of new hybrid schemes like FD-DG which would increase the accuracy in the unstructured part.
• Combining the SBP-SAT technique in space and time to obtain stable fully discrete approximations of IBVP.
We also mention the extension of SBP-SAT schemes to non-linear conservation laws (in the non-linear regime). The development of non-linear stability and convergence theory is of course hampered by the lack of mathematical knowledge. Existence results which are crucial when proving convergence are scarce for non-linear PDEs. However, the SBP-SAT schemes have many important properties that make them promising candidates in the developlment of non-linear stability theory and some steps have been taken in this direction.

References
[AC00] S. Abarbanel and A. E. Chertock. Strict stability of high-order compact implicit ﬁnite-diﬀerence schemes: The role of boundary conditions for hyperbolic PDEs, i. J. Comput. Phys., 160:42–66, 2000.
[ACY00] S. Abarbanel, A. E. Chertock, and A. Yefet. Strict stability of highorder compact implicit ﬁnite-diﬀerence schemes: The role of boundary conditions for hyperbolic PDEs, ii. J. Comput. Phys., 160:67–87, 2000.
[AD97] S. Abarbanel and A. Ditkowski. Asymptotically stable fourth-order accurate schemes for the diﬀusion equation on complex shapes. J. Comput. Phys., 133:279–288, 1997.
[AD99] S. Abarbanel and A. Ditkowski. Multi-dimensional asymptotically stable ﬁnite diﬀerence schemes for the advection-diﬀusion equation. Comput. & Fluids, 29:481–510, 1999.
[ADG00] S. Abarbanel, A. Ditkowski, and B. Gustafsson. On Error Bounds of Finite Diﬀerence Approximations to Partial Diﬀerential EquationsTemporal Behavior and Rate of Convergence. Journal of Scientiﬁc Computing, 2000.
[AG81] S. Abarbanel and D. Gottlieb. Optimal Time Splitting for Two- and Three-Dimensional Navier-Stokes Equations with Mixed Derivatives. J. Comput. Phys., 41:1–33, 1981.

22

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

[AN13] D. Amsallem and J. Nordstro¨m. High-order accurate diﬀerence schemes for the Hodgkin-Huxley equations. J. Comput. Phys., 252:573–590, 2013.
[BCVK02] H. Bijl, M.H. Carpenter, V.N. Vatsa, and C.A. Kennedy. Implicit time integration schemes for the unsteady compressible Navier-Stokes equations: laminar ﬂow. J. Comput. Phys., 179(1):313–329, 2002.
[BF88] F. Du Bois and P. Le Floch. Boundary conditions for nonlinear hyperbolic systems of conservation laws. J. Diﬀ. Eqs, 71(1):93–122, 1988.
[BN11] J. Berg and J. Nordstro¨m. Stable Robin solid wall boundary conditions for the Navier-Stokes equations. J. Comput. Phys., 230(19):7519–7532, 2011.
[BN12a] J. Berg and J. Nordstro¨m. Spectral analysis of the continuous and discretized heat and advection equation on single and multiple domains. Appl. Numer. Math., 62(11):1620–1638, 2012.
[BN12b] J. Berg and J. Nordstro¨m. Superconvergent functional output for time-dependent problems using ﬁnite diﬀerences on summation-byparts form. J. Comput. Phys., 231(20):6846–6860, 2012.
[BN13a] J. Berg and J. Nordstro¨m. Duality based boundary conditions and dual consistent ﬁnite diﬀerence discretizations of the Navier-Stokes and Euler equations. Technical report, under review in Journal of Computational Physics, 2013.
[BN13b] J. Berg and J. Nordstro¨m. On the impact of boundary conditions on dual consistent ﬁnite diﬀerence discretizations. J. Comput. Phys., 236(0):41–55, 2013.
[Bod09] D. J. Bodony. Scattering of an entropy disturbance into sound by a symmetric thin body. Physics of Fluids, 21, 2009.
[Bod10] D. J. Bodony. Accuracy of the Simultaneous-Approximation-Term boundary condition for time-dependent problems. J. Sci. Comput, 43(1):118–133, 2010.
[CGA94] M. H. Carpenter, D. Gottlieb, and S. Abarbanel. Time-stable boundary conditions for ﬁnite-diﬀerence schemes solving hyperbolic systems: Methodology and application to high-order compact schemes. J. Comput. Phys., 111(2), 1994.
[CKB+05] M. H. Carpenter, C.A. Kennedy, H. Bijl, S.A. Viken, and V. N. Vatsa. Fourth-order Runge-Kutta schemes for ﬂuid mechanics applications. J. Sci. Comput., 25:157–194, 2005.
[CNG99] M. H. Carpenter, J. Nordstro¨m, and D. Gottlieb. A Stable and Conservative Interface Treatment of Arbitrary Spatial Accuracy. J. Comput. Phys., 148, 1999.
[CNG10] M. H. Carpenter, J. Nordstro¨m, and D. Gottlieb. Revisiting and extending interface penalties for multi-domain summation-by-parts operators. J. Sci. Comput., 45(1-3):118–150, 2010.
[DBD+06] E. N. Dorband, E. Berti, P. Diener, E. Schnetter, and M. Tiglio. A numerical study of the quasinormal mode excitation of Kerr black holes. Phys. Rev. D, 75, 2006.
[DDS+07] P. Diener, E. N. Dorband, E. Schnetter, , and M. Tiglio. Optimized high-order derivative and dissipation operators satisfying summation by parts, and applications in three-dimensional multi-block evolutions. J. Sci. Comput., 32(1):109–145, 2007.
[EAN11] S. Eriksson, Q. Abbas, and J. Nordstro¨m. A stable and conservative method for locally adapting the design order of ﬁnite diﬀerence

SUMMATION-BY-PARTS

23

schemes. J. Comput. Phys., 230(11):4216–4231, 2011. [FC13] T. C. Fisher and M. H. Carpenter. High-order entropy stable ﬁnite
diﬀerence schemes for nonlinear conservation laws: Finite domains. J. Comput. Phys., 252:518–557, 2013. [FCN+13] T. C. Fisher, M. H. Carpenter, J. Nordstro¨m, N. K. Yamaleev, and C. Swanson. Discretely conservative ﬁnite-diﬀerence formulations for nonlinear conservation laws in split form: Theory and boundary conditions. J. Comput. Phys., 234(1):353–375, 2013. [FCYF11] T. C. Fisher, M. H. Carpenter, N. K. Yamaleev, and S.H. Frankel. Boundary closures for fourth-order energy stable weighted essentially non-oscillatory ﬁnite-diﬀerence schemes. J. Comput. Phys., 230:3727–3752, 2011. [FJAvW11] J. A. Fike, S. Jongsma, J.J. Alonso, and E. v.d. Weide. Optimization with gradient and hessian information calculated using hyper-dual numbers. In AIAA paper 2011-3807, 2011. [Gas13] G. Gassner. A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT ﬁnite diﬀerence methods. SIAM J. Sci. Comput., 35(3):A1233–A1253, 2013. [GKO95] B. Gustafsson, H.-O. Kreiss, and J. Oliger. Time dependent problems and diﬀerence methods. John Wiley & Sons, Inc., 1995. [GN07] J. Gong and J. Nordstro¨m. A stable and eﬃcient hybrid scheme for viscous problems in complex geometries. J. Comput. Phys., 226(2):1291–1309, 2007. [GN11] J. Gong and J. Nordstro¨m. Interface procedures for ﬁnite diﬀerence approximations of the advection-diﬀusion equation. J. Comput. Phys., 236:602–620, 2011. [GO96] M. Gerritsen and P. Olsson. Designing an eﬃcient solution strategy for ﬂuid ﬂows: 1. A stable high order ﬁnite diﬀerence scheme and sharp shock resolution for the Euler equations. J. Comput. Phys., 129(2):245–262, Dec. 1996. [GO98] M. Gerritsen and P. Olsson. Designing an eﬃcient solution strategy for ﬂuid ﬂows: II. Stable high-order central ﬁnite diﬀerence schemes on composite adaptive grids with sharp shock resolution. J. Comput. Phys., 147(2):293–317, Dec. 1998. [Gon06] J. Nordstro¨m & J. Gong. A stable hybrid method for hyperbolic problems. J. Comput. Phys., 212:436–453, 2006. [Gus75] B. Gustafsson. The convergence rate for diﬀerence approximations to mixed initial boundary value problems. Math. Comp., 29(130):396– 406, Apr. 1975. [Gus81] B. Gustafsson. The convergence rate for diﬀerence approximations to general mixed initial boundary value problems. SIAM J. Numer. Anal., 18(2):179–190, Apr. 1981. [Hic12] J.E. Hicken. Output error estimation for summation-by-parts ﬁnitediﬀerence schemes. J. Comput. Phys., pages 3828–3848, 2012. [HMIM07] F. Ham, K. Mattsson, G. Iaccarino, and P. Moin. Towards timestable and accurate les on unstructured grids. In Stavros C. Kassinos, Carlos A. Langer, Gianluca Iaccarino, and Parviz Moin, editors, Complex Eﬀects in Large Eddy Simulations, volume 56 of Lecture Notes in Computational Science and Engineering, pages 235–249. Springer Berlin Heidelberg, 2007. [HZ11] J.E. Hicken and D.W. Zingg. Superconvergent functional estimates from summation-by-parts ﬁnite-diﬀerence discretizations. SIAM

24

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

J.Sci.Comput., pages 893–922, 2011. [HZ13] J.E. Hicken and D.W. Zingg. Summation-by-parts operators and
high-order quadrature. J. Comput. Appl. Math, pages 111–125, 2013. [HZ14] J. E. Hicken and D. W. Zingg. Dual consistency and functional
accuracy: A ﬁnite-diﬀerence perspective. J. Comput. Phys., 256:161– 182, 2014. [JS96] G.-S. Jiang and C.-W. Shu. Eﬃcient implementation of weighted ENO schemes. J. Comput. Phys., 126:202–228, 1996. [KC03] C.A. Kennedy and M. H. Carpenter. Additive Runge-Kutta schemes for convection-diﬀusion-reaction equations. Appl. Numer. Math., 44:139–181, 2003. [KCL00] C.A. Kennedy, M. H. Carpenter, and R.M. Lewis. Low-storage, explicit Runge-Kutta schemes for the compressible Navier-Stokes equations. Appl. Numer. Math., 35, 2000. [KDN12] Jeremy E. Kozdon, Eric M. Dunham, and Jan Nordstro¨m. Interaction of Waves with Frictional Interfaces Using Summation-by-Parts Diﬀerence Operators: Weak Enforcement of Nonlinear Boundary Conditions. J. Sci. Comput., 50(2):341–367, 2012. [KDN13] J.E. Kozdon, E.M. Dunham., and J. Nordstro¨m. Simulation of dynamic earthquake ruptures in complex geometries using high-order ﬁnite diﬀerence methods. J. Sci. Comput., 55(1):92–124, 2013. [KL89] H.-O. Kreiss and J. Lorenz. Initial Boundary Value Problems and the Navier–Stokes Equations. Academic Press, New York, 1989. [KMRS09] U. Koley, S. Mishra, N.H. Risebro, and M. Sv¨ard. Higher order ﬁnite diﬀerence schemes for the magnetic induction equations. BIT, 49:375–395, 2009. [KMRS12] U. Koley, S. Mishra, N.H. Risebro, and M. Sv¨ard. Higher order ﬁnite diﬀerence schemes for the resistive magnetic induction equations. IMA Journal of Numerical Analysis, 32:1173–1193, 2012. [KO72] H.-O. Kreiss and J. Oliger. Comparison of accurate methods for the integration of hyperbolic equations. Tellus XXIV, 3, 1972. [KS74] H.-O. Kreiss and G. Scherer. Finite element and ﬁnite diﬀerence methods for hyperbolic partial diﬀerential equations. Mathematical Aspects of Finite Elements in Partial Diﬀerential Equations., Academic Press, Inc., 1974. [KS77] H.-O. Kreiss and G. Scherer. On the existence of energy estimates for diﬀerence approximations for hyperbolic systems. Technical report, Dept. of Scientiﬁc Computing, Uppsala University, 1977. [KW93] H.-O. Kreiss and L. Wu. On the stability deﬁnition of diﬀerence approximations for the initial boundary value problem. Applied Numerical Mathematics, 12:213–227, 1993. [LN10] J. Lindstr¨om and J. Nordstro¨m. A stable and high-order accurate conjugate heat transfer problem. J. Comput. Phys., 229(14):5440– 5456, 2010. [LN13] T. Lundquist and J. Nordstro¨m. The SBP-SAT technique for timediscretization. In AIAA paper 2013-2834, 21st AIAA Computational Fluid Dynamics Conference, 2013. [LR56] P.D. Lax and R.D. Richtmyer. Survey of the stability of linear ﬁnite diﬀerence equations. Comm. on Pure and Applied Math., IX, 1956. [LT98] D. Levy and E. Tadmor. From semi-discrete to fully-discrete: stability of Runge-Kutta schemes by the energy method. SIAM Review, 40:40–73, 1998.

SUMMATION-BY-PARTS

25

[Mat03] K. Mattsson. Boundary procedures for summation-by-parts operators. J. Sci. Comput., 18, 2003.
[Mat11] K. Mattsson. Summation by parts operators for ﬁnite diﬀerence approximations of second-derivatives with variable coeﬃcients. Journal of Scientiﬁc Computing, pages 1–33, 2011.
[MC10] K. Mattsson and M. H. Carpenter. Stable and accurate interpolation operators for high-order multiblock ﬁnite diﬀerence methods. SIAM J. Sci. Comput., 32(4):2298–2320, 2010.
[MHI08] K. Mattsson, F. Ham, and G. Iaccarino. Stable and accurate wavepropagation in discontinuous media. J. Comput. Phys., 227:249–269, 2008.
[MHI09] K. Mattsson, F. Ham, and G. Iaccarino. Stable boundary treatment for the wave equation on second-order form. J. Sci. Comput., 41(3):366–383, 2009.
[Mic83] D. Michelson. Stability theory of diﬀerence approximations for multidimensional initial-boundary value problems. Math. Comp., 40:1–45, 1983.
[Mic87] D. Michelson. Convergence theorem for diﬀerence approximations of hyperbolic quasi-linear initial-boundary value problems. Math. Comp., 49:445–459, 1987.
[MMR11] M. Motamed, C. B. Macdonald, and S. J. Ruuth. On the linear stability of the ﬁfth-order WENO discretization. J. Sci. Comput., 47(2):127–149, 2011.
[MN04] K. Mattsson and J. Nordstro¨m. Summation by parts operators for ﬁnite diﬀerence approximations of second derivatives. J. Comput. Phys., 199(2), 2004.
[MN06] K. Mattsson and J. Nordstro¨m. High order ﬁnite diﬀerence methods for wave propagation in discontinuous media. J. Comput. Phys., 220:249–269, 2006.
[MP10] K. Mattsson and F. Parisi. Stable and accurate second-order formulation of the shifted wave equation. Comm. Comp. Phys., 7:101–137, 2010.
[MS10] S. Mishra and M. Sv¨ard. On stability of numerical schemes via frozen coeﬃcients and the magnetic induction equations. BIT, 50:85–108, 2010.
[MSCN07] K. Mattsson, M. Sv¨ard, M.H. Carpenter, and J. Nordstro¨m. Highorder accurate computations for unsteady aerodynamics. Comput. & Fluids, 36(3):636–649, 2007.
[MSN04] K. Mattsson, M. Sv¨ard, and J. Nordstro¨m. Stable and Accurate Artiﬁcial Dissipation. J. Sci. Comput., 21(1):57–79, August 2004.
[MSS08] K. Mattsson, M. Sv¨ard, and M. Shoeybi. Stable and accurate schemes for the compressible Navier-Stokes equations. J. Comput. Phys., 227:2293–2316, 2008.
[Mu¨l08] B. Mu¨ller. High order numerical simulation of aeolian tones. Comput. & Fluids, 37(4):450–461, 2008.
[MY02] B. Mu¨ller and H. Yee. Entropy splitting for high order numerical simulation of vortex sound at low Mach numbers. J. Sci. Comput., 17:181–190, 2002.
[MY09] B. Mu¨ller and H. Yee. Numerical simulation of conﬁned pulsating jets in human phonation. Comput. & Fluids, 38(7):1375–1383, 2009.
[NB01] J. Nordstro¨m and M. Bj¨orck. Finite volume approximations and strict stability for hyperbolic problems. Appl. Numer. Math.,

26

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

38(3):237–255, 2001. [NB13] J. Nordstro¨m and J. Berg. Conjugate heat transfer for the unsteady
compressible Navier-Stokes equations using a multi-block coupling. Comput. & Fluids, 72:20–29, 2013. [NC99] J. Nordstro¨m and M. H. Carpenter. Boundary and interface conditions for high-order ﬁnite-diﬀerence methods applied to the Euler and Navier-Stokes equations. J. Comput. Phys., 148, 1999. [NC01] J. Nordstro¨m and M. H. Carpenter. High-order ﬁnite diﬀerence methods, multidimensional linear problems, and curvilinear coordinates. J. Comput. Phys., 173, 2001. [NE10] J. Nordstro¨m and S. Eriksson. Fluid structure interaction problems: The necessity of a well posed, stable and accurate formulation. Communications in Computational Physics, 8(5):1111–1138, 2010. [NEE12] J. Nordstro¨m, S. Eriksson, and P. Eliasson. Weak and Strong Wall Boundary Procedures and Convergence to Steady-State of the Navier-Stokes Equations. J. Comput. Phys., 231(14):4867–4884, 2012. [NFAE03] J. Nordstro¨m, K. Forsberg, C. Adamsson, and P. Eliasson. Finite volume methods, unstructured meshes and strict stability for hyperbolic problems. Appl. Numer. Math., 45(4):453–473, June 2003. [NG03] J. Nordstro¨m and R. Gustafsson. High order ﬁnite diﬀerence approximations of electromagnetic wave propagation close to material discontinuities. J. Sci. Compt., 18(2):215–234, 2003. [NGvdWS09] J. Nordstro¨m, J. Gong, E. van der Weide, and M. Sv¨ard. A stable and conservative high order multi-block method for the compressible Navier-Stokes equations. J. Comput. Phys., 228:9020–9035, 2009. [NHS+09] J. Nordstro¨m, F. Ham, M. Shoeybi, E. van der Weide, M. Sv¨ard, K. Mattsson, G. Iaccarino, and J. Gong. A hybrid method for unsteady inviscid ﬂuid ﬂow. Comput. & Fluids, 38:875–882, 2009. [NKG12] A. Nissen, G. Kreiss, and M. Gerritsen. Stability at nonconforming grid interfaces for a high order discretization of the schr¨odinger equation. Journal of Scientiﬁc Computing, 53(3):528–551, 2012. [NL13] J. Nordstro¨m and T. Lundquist. Summation-by-parts in time. J. Comput. Phys., 251:487–499, 2013. [NMS07] J. Nordstro¨m, K. Mattsson, and C. Swanson. Boundary conditions for a divergence free velocity-pressure formulation of the NavierStokes equations. J. Comput. Phys., 225(1):874–890, 2007. [Nor95] J. Nordstro¨m. The Use of Characteristic Boundary Conditions for the Navier-Stokes Equations. Computers & Fluids, 24(5):609–623, 1995. [Nor06] J. Nordstro¨m. Conservative ﬁnite diﬀerence formulation, variable coeﬃcients, energy estimates and artiﬁcial dissipation. J. Sci. Comput., 29(3):375–404, 2006. [Nor07] J. Nordstro¨m. Error bounded schemes for time-dependent hyperbolic problems. SIAM J. Sci. Comput., 30(1):46–59, 2007. [Ols95a] P. Olsson. Summation by parts, projections, and stability I. Math. Comp., 64:1035, 1995. [Ols95b] P. Olsson. Summation by parts, projections, and stability II. Math. Comp., 64:1473, 1995. [OO94] P. Olsson and J. Oliger. Energy and maximum norm estimates for nonlinear conservation laws. Technical report 94.01, The Research Institute of Advanced Computer Science, Jan. 1994.

SUMMATION-BY-PARTS

27

[PC81] T.H. Pulliam and D.S. Chaussee. A diagonal form of an implicit approximate-factorization algorithm. J. Comput. Phys., 39:347–363, 1981.
[PDN13] P. Pettersson, A. Doostan, and J. Nordstro¨m. On stability and monotonicity requirements of ﬁnite diﬀerence approximations of stochastic conservation laws with random viscosity. Computer Methods in Applied Mechanics and Engineering, 258:134–151, 2013.
[PIN09] P. Pettersson, G. Iaccarino, and J. Nordstro¨m. Numerical analysis of the Burgers’ equation in the presence of uncertainty. J. Comput. Phys., 228(22):8394–8412, 2009.
[PIN13] P. Pettersson, G. Iaccarino, and J. Nordstro¨m. An intrusive hybrid method for discontinuous two-phase ﬂow under uncertainty. Comput. & Fluids, 86:228–239, 2013.
[PIN14] P. Pettersson, G. Iaccarino, and J. Nordstro¨m. A stochastic Galerkin method for the Euler equations with Roe variable transformation. Journal of Computational Physics, 257:481–500, 2014.
[PNI10] P. Pettersson, J. Nordstro¨m, and G. Iaccarino. Boundary procedures for the time-dependent Burgers’ equation under uncertainty. Acta Mathematica Scientia, 30(2):539–550, 2010.
[SCN07] M. Sv¨ard, M.H. Carpenter, and J. Nordstro¨m. A stable high-order ﬁnite diﬀerence scheme for the compressible Navier-Stokes equations, far-ﬁeld boundary conditions. J. Comput. Phys., 225:1020–1038, 2007.
[SDDT06] E. Schnetter, P. Diener, E. N. Dorband, and M. Tiglio. A multiblock infrastructure for three-dimensional time-dependent numerical relativity. Class. Quantum Grav, 23:S553–S578, 2006.
[SGN06] M. Sv¨ard, J. Gong, and J. Nordstro¨m. Stable artiﬁcial dissipation operators for ﬁnite volume schemes on unstructured grids. Appl. Numer. Math., 56(12):1481–1490, Dec. 2006.
[SLN10] M. Sv¨ard, J. Lundberg, and J. Nordstro¨m. A computational study of vortex-airfoil interaction using high-order ﬁnite diﬀerence methods. Comput. & Fluids, 39:1267–1274, 2010.
[SM11] M. Sv¨ard and S. Mishra. Implicit-explicit schemes for ﬂow equations with stiﬀ source terms. Journal of Computational and Applied Mathematics, 235:1564–1577, 2011.
[SM12] M. Sv¨ard and S. Mishra. Entropy stable schemes for initial-boundaryvalue conservation laws. ZAMP, 63:985–1003, 2012.
[SMN05] M. Sv¨ard, K. Mattsson, and J. Nordstro¨m. Steady state computations using summation-by-parts operators. J. Sci. Comput., 24(1):79–95, September 2005.
[SN04] M. Sv¨ard and J. Nordstro¨m. Stability of ﬁnite volume approximations for the Laplacian operator on quadrilateral and triangular grids. Appl. Numer. Math., 51(1):101–125, October 2004.
[SN06] M. Sv¨ard and J. Nordstro¨m. On the order of accuracy for diﬀerence approximations of initial-boundary value problems. J. Comput. Physics, 218(1):333–352, October 2006.
[SN08] M. Sv¨ard and J. Nordstro¨m. A stable high-order ﬁnite diﬀerence scheme for the compressible Navier-Stokes equations, no-slip wall boundary conditions. J. Comput. Phys., 227:4805–4824, 2008.
[SO¨ 13] M. Sv¨ard and H. O¨ zcan. Entropy stable schemes for the Euler equations with far-ﬁeld and wall boundary conditions. online ﬁrst, Journal of Scientiﬁc Computing, 2013.

28

MAGNUS SVA¨ RD AND JAN NORDSTRO¨ M

[SSHM10] M. Shoeybi, M. Sv¨ard, F. Ham, and P. Moin. An adaptive implicitexplicit scheme for the DNS and LES of compressible ﬂows on unstructured grids. J. Comput. Phys., 229:5944–5965, 2010.
[Str64] G. Strang. Accurate partial diﬀerence methods II. non-linear problems. Num. Math., 6:37–46, 1964.
[Str94] B. Strand. Summation by Parts for Finite Diﬀerence Approximations for d/d x. J. Comput. Phys., 110, 1994.
[Sv¨a04] M. Sv¨ard. On coordinate transformation for summation-by-parts operators. J. Sci. Comput., 20(1), 2004.
[Sv¨a12] M. Sv¨ard. Third-order accurate entropy-stable scheme for initialboundary-value conservation laws. ZAMP, 63:599–623, 2012.
[Sv¨a13] M. Sv¨ard. A note on l∞ bounds and convergence rates of summationby-parts schemes. Technical report, under review, BIT, 2013.
[SW88] J. C. Strikwerda and B. A. Wade. An extension of the Kreiss matrix theorem. SIAM J. Numer. Anal., 25(6):1272–1278, 1988.
[Tad03] E. Tadmor. Entropy stability theory for diﬀerence approximations of nonlinear conservation laws and related time-dependent problems. Acta Numerica, pages 451–512, 2003.
[vdWS12] E. van der Weide and M. Sv¨ard. Test cases c1.1, c1.2, c1.3. In 1st International Workshop on High-Order CFD Methods January 7-8, 2012. Nashville, Tennessee, USA, 2012. http://zjw.public.iastate.edu/hiocfd/abstracts/C11 Weide.pdf / http://zjw.public.iastate.edu/hiocfd/summary/C1.1 summary.pdf / http://zjw.public.iastate.edu/hiocfd/summary/C1.2 summary.pdf / http://zjw.public.iastate.edu/hiocfd/summary/C1.3 summary.pdf.
[Wad90] B. A. Wade. Symmetrizable ﬁnite diﬀerence operators. Math. Comp, 54:525–543, Apr. 1990.
[WCVJ13] D. M. Williams, P. Castonguay, P. E. Vincent, and A. Jameson. Energy stable ﬂux reconstruction schemes for advection-diﬀusion problems on triangles. J. Comput. Phys., 250:53–76, 2013.
[WS07] R. Wang and R. J. Spiteri. Linear instability of the ﬁfth-order WENO method. SIAM J. Numer. Anal., 45(5):1871–1901, 2007.
[YC09] N.K. Yamaleev and M.H. Carpenter. Third-order energy stable WENO scheme. J. Comput. Phys., 228(8):3025–3047, 2009.
[ZST08] B. Zink, E. Schnetter, and M. Tiglio. Multi-patch methods in general relativistic astrophysics - i. hydrodynamical ﬂows on ﬁxed backgrounds. Phys. Rev. D, 2008, 2008.

