Optimal diagonal-norm SBP operators

Journal of Computational Physics 264 (2014) 91–111

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Optimal diagonal-norm SBP operators Ken Mattsson 1 , Martin Almquist, Mark H. Carpenter Department of Information Technology, Uppsala University, PO Box 337, S-751 05 Uppsala, Sweden

a r t i c l e

i n f o

Article history: Received 28 June 2013 Received in revised form 13 December 2013 Accepted 19 December 2013 Available online 15 January 2014 Keywords: Finite difference methods Multi-block grids High order accuracy Stability Boundary treatment Euler equations

a b s t r a c t Optimal boundary closures are derived for first derivative, finite difference operators of order 2, 4, 6 and 8. The closures are based on a diagonal-norm summation-by-parts (SBP) framework, thereby guaranteeing linear stability on piecewise curvilinear multi-block grids and entropy stability for nonlinear equations that support a convex extension. The new closures are developed by enriching conventional approaches with additional boundary closure stencils and non-equidistant grid distributions at the domain boundaries. Greatly improved accuracy is achieved near the boundaries, as compared with traditional diagonalnorm operators of the same order. The superior accuracy of the new optimal diagonal-norm SBP operators is demonstrated for linear hyperbolic systems in one dimension and for the nonlinear compressible Euler equations in two dimensions. © 2014 Elsevier Inc. All rights reserved.

1. Introduction Wave-propagation problems frequently require farfield boundaries to be positioned many wavelengths away from the disturbance source. Efficient simulation of these problems requires numerical techniques capable of accurately propagating disturbances over long distances. It is well known that high-order finite difference methods (HOFDMs) are ideally suited for problems of this type. (See the pioneering paper by Kreiss and Oliger [24].) Not all high-order spatial operators are applicable, however. For example, schemes that are G-K-S stable [15], while being convergent to the true solution as x → 0, may experience nonphysical solution growth in time [5], thereby limiting their efficiency for long-time simulations. Thus, it is imperative to use HOFDMs that do not allow nonphysical solution growth in time—a property termed “strict stability” [14]. Deriving a strictly stable, accurate, and conservative HOFDM is a significant challenge that has received considerable past attention. (For examples, see [26,44,41,1,3,17,42,16].) A robust and well-proven high-order finite difference methodology that ensures the strict stability of time-dependent partial differential equations (PDEs) is the summation-by-parts–simultaneous approximation term (SBP-SAT) method. The SBP-SAT method combines semi-discrete operators that satisfy a summation-by-parts (SBP) formula [22], with physical boundary conditions implemented using the simultaneous approximation term (SAT) method [5]. (The SAT method was originally developed for pseudo-spectral approximations [11,12]. Examples of the pseudo-spectral penalty approach can be found in [17,50].) Examples of the SBP-SAT approach can be found in [36–38,30,32,33,39,29,46,25,9,31,19,18,21]. An added benefit of the SBP-SAT method is that it naturally extends to multi-block geometries while retaining the essential single-block properties: strict stability, accuracy, and conservation [6]. Thus, problems involving complex domains or non-smooth geometries are easily amenable to the approach. Refs. [29,32,20] report applications of the SBP-SAT method to problems involving nontrivial geometries.

1

E-mail address: [email protected] (K. Mattsson). Tel.: +46 18 4717631; fax: +46 18 523049.

0021-9991/$ – see front matter © 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.12.041

92

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

A wide variety of numerical methods may be expressed as summation-by-parts operators; examples include finite difference operators, spectral collocation [8,7], and some finite volume methods [48]. Herein we focus exclusively on finite difference operators. Finite difference operators may be further categorized by the structure of their norm: a) diagonal, b) diagonal interior with block boundary closures, c) fully banded. The early papers [23,43] report first-derivative SBP operators up to 6th order, based on diagonal-norms and block-norms, while in [4] a 4th-order accurate first derivative Padé type banded SBP operator is derived. The SBP framework may also be used to develop second-derivative operators. In [30] diagonal- and block-norm SBP operators up to 8th order were derived for both the first and second derivatives. Matrix structural constraints required by all SBP operators preclude uniform spatial accuracy for diagonal-norm finite difference operators [23]. Indeed, a 2pth order diagonal-norm SBP operator may be closed at the boundaries with at most pth-order accurate biased stencils. For first order hyperbolic problems, this implies a suboptimal convergence rate (i.e., global convergence) of ( p + 1)th order when using a 2pth order interior operator. In contrast, a block-norm SBP operator with an identical interior stencil may be closed at the boundaries with (2p − 1)th-order accurate stencils, thereby yielding optimal 2pth order convergence. Similarly, strongly parabolic problems and second order hyperbolic problems converge at a rate of ( p + 2)th order using a 2pth order diagonal-norm SBP operator. (See [13,49] for more information on the accuracy of finite difference approximations.) At first glance, it would appear that block-norm or banded SBP operators would be preferred for simulations with stringent accuracy requirements. It is well known, however, that block-norm SBP operators may exhibit instabilities when applied to problems mapped onto curvilinear grids [45], while diagonal-norm operators remain stable. The same instability mechanism exists on Cartesian grids for second order systems with variable coefficients [34]. Furthermore, recent developments appearing in Ref. [10] provide a general procedure for developing entropy-consistent and entropy-stable SBP schemes for systems of equations that support a convex extension (e.g., the entropy in the compressible Navier–Stokes equations). This procedure is only valid, however, for diagonal-norm SBP spatial operators. For these reasons diagonal-norm operators have gained more interest in practical applications. State-of-the-art diagonal-norm finite-difference operators approach a point of diminishing return with 8th order interior operators and 4th order boundary closures [30,32], thereby achieving 5th order convergence (for first order hyperbolic systems and incompletely parabolic systems) on curvilinear grids. A 10th order diagonal-norm SBP operator was constructed in [9] by extending the number of boundary stencils, but exhibits large negative eigenvalues that separate from the cluster, leading to very small time steps for explicit time-stepping methods. Herein, diagonal-norm SBP finite difference operators are derived with highly optimized boundary closures. Although the resulting operators are formally 5th-order accurate (for the 8th order case), the leading order error constants are reduced by several orders of magnitude. Thus, very accurate simulations can be performed on relatively coarse meshes. The optimized boundary operators utilize additional unknowns in the SBP boundary closures, as well as a non-equidistant point distribution at the first three boundary points in the domain. The superior accuracy properties of the newly derived diagonal-norm SBP operators are demonstrated for linear hyperbolic problems in 1-D and the compressible Euler equations in 2-D, where the computational domain consists of two curvilinear grids. In Section 2 the SBP-SAT method is introduced in one spatial dimension. In Section 3 the accuracy and stability properties of the newly developed SBP operators are verified by performing 1-D numerical simulations. Verification of accuracy and stability by numerical studies of the 2-D problem is performed in Section 4. Section 5 summarizes the work. The SBP operators are presented in Appendices A and B. 2. The finite difference method The focus in the present study is the derivation of accuracy-optimized diagonal-norm SBP operators. The traditional diagonal-norm SBP operators first presented in [22,23] and later re-derived using Maple in [43,30] suffer from low accuracy at the boundaries. As a remedy one can introduce a non-diagonal norm, which allows for much more accurate boundary closures. However, non-diagonal SBP operators suffer from stability issues related to curvilinear grids [45] and variablecoefficient problems [35]. Furthermore, the recent development in [10]—concerning entropy-stable SBP schemes for the compressible Navier–Stokes equations—requires diagonal-norm SBP operators. 2.1. Definitions The following definitions are needed later in the present study. Let u, v ∈ L 2 [0, 1] where u T = [u (1) , u (2) , . . . , u (k) ] and v = [ v (1) , v (2) , . . . , v (k) ] are real-valued vector functions with k components. Let the inner product be defined by (u, v) = 1 T u A(x)v dx, A(x) = A T (x) > 0, and let the corresponding norm be u2A = (u, u). 0 The domain (0  x  1) is discretized using m grid points: T

x = [x1 , x2 , . . . , xm−1 , xm ] T , i.e., x denotes a vector holding the grid-points. For an equidistant grid we have

xi = (i − 1)h,

i = 1, 2, . . . , m , h =

1 m−1

.

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

93

In the present study we will allow the 3 grid-points closest to the boundary to be non-equidistant to improve accuracy. ( j) ( j) ( j) The approximate solution vector is given by v T = [ v (1) , v (2) , . . . , v (k) ], where v ( j ) = [ v 1 , v 2 , . . . , v m ] is the discrete solution vector of the jth component. Similarly, we define an inner product for discrete real-valued vector functions u , v ∈ ˜ v, where H˜ = I k ⊗ H is positive definite (H is a symmetric positive definite m × m matrix) and Rk×m by (u , v ) H A = u T A H A is the projection of A(x) onto the block diagonals. If A(x) is a diagonal k × k matrix, A is a diagonal (km) × (km) matrix. ˜ A v. I k denotes the unit matrix of size k × k. The corresponding norm is  v 2H A = v T H

˜ A defines a norm if and only if H˜ A is positive definite. For variable A (x) this can only be Remark. The matrix product H guaranteed if the norm H is a diagonal matrix (see [45] for a detailed study on this). Variable coefficients are often present due to the underlying physics or the geometry requiring the introduction of curvilinear grids. We will focus on both these issues in the present study. Here we make use of the Kronecker product:

⎡ ⎢

C ⊗D =⎣

c 0,0 D

.. . c p −1,0 D

···

c 0,q−1 D

.. .

⎤

⎥ ⎦,

· · · c p −1,q−1 D

where C is a p × q matrix and D is an m × n matrix. Two useful rules for the Kronecker product are ( A ⊗ B )(C ⊗ D ) = ( AC ) ⊗ ( B D ) and ( A ⊗ B ) T = A T ⊗ B T . The following vectors will be frequently used:

e 1 = [1, 0, . . . , 0] T ,

em = [0, . . . , 0, 1] T .

(1)

The following definition (first stated in [33,29]) is central to the present study: Definition 2.1. An explicit 2pth-order accurate finite difference scheme with minimal stencil width for the Cauchy problem is denoted a 2pth-order accurate narrow-stencil. 2.2. Traditional SBP operators The SBP definitions can be found in earlier papers (see for example [30,33,29,28,31,34]). For completeness we restate them below for the first derivative: T Definition 2.2. A difference operator D 1 = H −1 ( Q − 12 e 1 e 1T + 12 em em ), approximating ∂/∂ x, using a 2pth-order accurate narrow-stencil, is said to be a 2pth-order first-derivative SBP operator if H is symmetric and positive definite, and Q + Q T = 0. T Definition 2.3. A difference operator D 1 = H −1 ( Q − 12 e 1 e 1T + 12 em em ), approximating ∂/∂ x, using a 2pth-order accurate narrow-stencil, is said to be a 2pth-order diagonal-norm first-derivative SBP operator if H is diagonal and positive definite, and Q + Q T = 0.

The following two new definitions are central to the present study. Definition 2.4. Let xq be the projection of the polynomial

xq q!

onto the discrete grid-points, i.e., a vector denoted x.

xq q!

onto the discrete grid-points, i.e., a vector denoted x.

T Let e (q) = Hxq−1 − ( Q − 12 e 1 e 1T + 12 em em )xq be the qth order error vector. We say that the diagonal-norm first-derivative 1 1 T SBP operator D 1 = H −1 ( Q − 2 e 1 e 1T + 2 em em ) is 2pth-order if e (q) vanishes for q = 1 . . . 2p, in the interior and for q = 1 . . . p, at the boundaries.

Definition 2.5. Let xq be the projection of the polynomial

T Let e (q) = Hxq−1 − ( Q − 12 e 1 e 1T + 12 em em )xq be the qth order error vector. The discrete l2 -norm of the error e (q) is defined 2 T as e (q) h = he ( p ) e ( p ) .

Remark. The boundary closure for a 2pth-order diagonal-norm first-derivative SBP operator is restricted to pth-order accuracy (see [30]) at the boundaries, i.e., e ( p +1) is the leading order error. This implies that the convergence rate drops to ( p + 1)th order for first order hyperbolic problems, and ( p + 2)th order for parabolic problems and second order hyperbolic problems. (See [13,49] for more information on the accuracy of finite difference approximations.) Constructing 2pth-order diagonal-norm first-derivative SBP operators requires a minimum of 2p non-central difference stencils (that are pth order accurate) at the boundaries. In this paper the diagonal-norm first-derivative SBP operators with

94

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

a minimum number of boundary stencils will be referred to as traditional, to distinguish them from the accuracy-optimized diagonal-norm SBP operators derived in the present study. Traditional first-derivative SBP operators were first constructed in [22,23] and later rederived using Maple in [43,30]. It can be shown [23] that traditional diagonal-norm first-derivative SBP operators do not exist beyond 8th order. The traditional 6th order diagonal-norm SBP operator has one “free” parameter in Q and the corresponding 8th order operator has 3 “free” parameters in Q . In [43] the free parameters were chosen to yield an SBP operator with minimum band-width at the boundaries (the motivation was based on implementation issues, not considering numerical efficiency). However, optimizing for band-width yields operators with poor numerical behaviour, which was first found in [32]. Numerical efficiency can mean many things, but here we will focus on two relevant measures: 1) the spectral radius of the operator, and 2) the leading order truncation error of the boundary stencils. A careful numerical study was performed in [32], minimizing the spectral radius (by tuning the free parameters) of the traditional 6th-order and 8th-order diagonal-norm SBP operators, which led to highly improved numerical efficiency. A strong correlation between optimal spectral radius and optimal leading order truncation error for the SBP operators was found in [33]. Hence, by minimizing the norm of the leading order truncation error e ( p +1) h2 (see Definition 2.5), we expect to also minimize the spectral radius of the SBP operator. (In [33] we further optimized for yet another important property, referred to as compatibility between first- and second-derivative SBP operators which are important in the presence of mixed second-derivative terms.) The traditional first-derivative SBP operators used in the present study have the optimal parameter choice concerning both accuracy and spectral radius, i.e., they have optimal numerical properties. This is achieved by choosing the “free” parameters in Q to minimize the discrete l2 -norm of the leading order error, i.e., minimizing e ( p +1) h2 for the 2pth order case (see Definition 2.5). Remark. The traditional 2pth order block-norm SBP operators have more free parameters than the corresponding diagonalnorm SBP operators of the same order, and the leading order error is given by e (2p −1) , resulting in full order convergence. In a recent study [27] up to 10th order traditional block-norm SBP operators were constructed, where the free parameters were chosen to minimize the leading order truncation error (and hence the spectral radius). 2.3. Extended diagonal-norm SBP operators Despite choosing the “free” parameters in the traditional diagonal-norm operators optimally, the accuracy of the boundary stencils remains poor compared to the internal central finite difference scheme. For the traditional 4th-order case this is not a major issue, but for the 6th-order case and the 8th-order case in particular, there is a large discrepancy between the internal truncation error and the leading order truncation error at the boundaries. In [9] a 10th order first-derivative diagonal-norm SBP operator was constructed by introducing 11 boundary stencils, resulting in 10 “free” parameters, one of which appearing in the norm H rather than Q . Here one needs choose the free parameters extra carefully; the SBP requirement that H should be positive definite must not be violated. In [9], the attempt to achieve an efficient 10th order diagonal SBP operator by minimizing the leading order truncation error failed, since a straightforward attempt will violate the necessary SBP condition H > 0. The resulting 10th order SBP operator was in fact less accurate than the corresponding 8th order operator (for realistic grid-refinements). In [27] a new attempt to build a 10th order diagonal-norm SBP operator with 11 boundary stencils was performed, employing a numerical minimization for both the leading order truncation error and the spectral radius, which resulted in improved accuracy compared to the 8th-order SBP operator, although the difference is small (see for example Table 4). It can be shown that any attempt to construct a diagonal-norm SBP operator with more than 8 boundary stencils leads to “free” parameters in H . The necessary condition H > 0 then severely restricts the parameter space and limits the optimization procedure. Thus, it seems difficult to construct efficient diagonal-norm SBP operators with more than 8 boundary stencils. To reduce the leading order truncation error, one can introduce more than the traditional 2p boundary stencils, which leads to more “free” parameters. This approach will only work up to the 6th-order case, since a valid SBP solution to this minimization cannot be obtained using more than 8 boundary points. In the present study we have constructed novel 2nd, 4th and 6th order diagonal-norm SBP operators with 2 additional boundary stencils (at each boundary), compared to the corresponding traditional SBP operators. The free parameters are chosen to minimize the discrete l2 -norm of the leading order truncation errors (see Definition 2.5). For the 6th-order case we first minimize e (4) h2 , which leaves one free parameter. The last parameter is determined by minimizing e (5) h2 . These newly derived SBP operators will be referred to as extended and are listed in Appendix A.

Remark. One might hope to construct an extended 6th-order operator that makes e (4) h2 vanish. Unfortunately this cannot be achieved. 2.4. Optimal diagonal-norm SBP operators A novel idea to further reduce the leading order truncation errors—also for the 8th-order diagonal-norm SBP operator—is to allow a non-equidistant grid at the boundaries. Hence, we allow the location of the first boundary points to be “free”

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

95

Fig. 1. The locations of the grid-points for the optimal 4th-, 6th- and 8th-order accurate cases, using m = 11. The locations of the first three grid-points (closest to the boundary) are non-equidistant. A Cartesian grid would have the grid-points at the vertical dotted lines. Table 1 The normalised spectral radius of the various SBP operators (normalised against m, to be independent of h). In comparison, the traditional 8th order SBP operator optimized for bandwidth yields a value of 15.71 (instead of 1.69). The 10th order SBP operator (presented in [27]) has a normalised spectral radius of 1.79. Operator

4th

6th

8th

Traditional Extended Optimal

1.34 1.34 1.36

1.55 1.55 1.52

1.69 1.63

in the minimization. Allowing for a “free” location of the boundary points, however, leads to non-linear equations when minimizing the discrete l2 -norm of the leading order truncation errors (see Definition 2.5). In the present study we make use of the numerical minimization routine Minimize in Maple. We have found that no more than three “free” boundary points can be allowed using the numerical minimization routine in Maple (it simply becomes too computationally demanding). It is possible that a more efficient minimization routine could allow for more “free” boundary points, but finding an optimal nonlinear minimization routine is out of scope in the present study. Let x ∈ [0, 1]. Introduce d = d1 + d2 + d3 , where hd1 is the distance between the first and second grid-point, hd2 is the distance between the second and third grid-point, and hd3 is the distance between the third and fourth grid-point. Here h = 1/(2d + (m − 7)) and m is the number of grid-points. The grid-points are located at









x = 0, d1 h, (d1 + d2 )h, dh, (d + 1)h, (d + 2)h, . . . , d + (m − 8) h, 1 − dh, 1 − (d1 + d2 )h, 1 − d1 h, 1 , where d1 , d2 and d3 are unknowns in the minimization. We have derived 2nd, 4th, 6th and 8th order accurate diagonal SBP operators, where the location of the first three grid-points (closest to the boundary) are non-equidistant, see Fig. 1. For the 2nd-order case we use three boundary stencils, and for the 4th and 6th-order cases we use 5 and 7 boundary stencils, respectively. For the 8th-order case we are restricted to 8 boundary stencils. In the present study these new diagonal-norm SBP operators will be referred to as optimal. The nonlinear equations are solved using a numerical minimization routine in Maple. There are many possible local minima, and we cannot guarantee global minimum. To limit the search space we have restricted the values of d1 , d2 and d3 to be between 0.3 and 1.5, to prevent small grid-spacings that could result in a large spectral radius of the resulting SBP operator. A careful numerical study of the resulting SBP operators indicate that it is favourable to minimize a linear combination of the two leading order errors, i.e., to minimize e ( p +1) + e ( p +2) h2 (instead of e ( p +1) h2 ). We also tried minimizing e ( p +1) +

e ( p +2) + e ( p +3) h2 , but that did not improve the numerical efficiency of the resulting SBP operator. In Table 1 we have listed the normalised spectral radii of the traditional, extended and optimal 4th, 6th and 8th order SBP operators. We have normalised against the number of grid-points m to yield a measure that is independent of h. For a given order, the differences between the operators are small. It is worth mentioning that the traditional 8th-order SBP operator has a normalised spectral radius of 1.69 when optimizing for the leading order truncation error, while it grows to 15.71 when optimized for bandwidth, as presented in [43]. All operators are presented in Appendices A and B. The results presented in Sections 3 and 4 show that there is a huge gain in accuracy by allowing the grid-points to be non-equidistant at the boundaries.

96

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

Remark. The process of deriving optimal SBP operators yields the optimal positions for the grid-points, resulting in a non-equidistant grid close to the boundaries. The gain in efficiency is rather high, especially for the 8th-order case. 2.5. The SBP-SAT method Consider

ut = Aux ,

0  x  1 , t  0,

where A = A T is a constant coefficient matrix. (We also assume initial data u(x, 0) = f(x).) This equation is now mapped onto a curvilinear grid x = x(ξ ),

Jut = Auξ ,

0  ξ  1, t  0,

(2)

where J = ξx−1 is the Jacobian matrix. We diagonalize A by a similar transformation, A = T T ΛT, where Λ holds the eigenvalues to A. We will apply boundary conditions to the ingoing characteristics, to specify the correct number of boundary conditions. We introduce





2A± = T T Λ ± |Λ| T.

(3)

At the left (L) and right (R) boundaries we specify data for A+ u and A− u respectively, i.e.,

A− u L = A− g L ,

A+ u R = A+ g R ,

(4)

where u L , R are the unknowns, and g L , R the corresponding data at the Left and Right boundaries, respectively. Multiplying (2) by u T , integrating by parts, adding the transpose and imposing the ingoing characteristics (4) leads to,

d dt

u2J = uTL A− u L + gTL A+ g L − uTR A+ u R − gTR A− g R .

(5)

Hence, the problem defined by (2) and (4) is strongly well posed [14]. The semi-discrete problem using the SBP-SAT method reads

˜J v t = A ⊗ D 1 + A − ⊗ H −1 e 1 ( v L − g L ) − A + ⊗ H −1 em ( v R − g R ),

(6)

where ˜J = I k × J and J is the projection of the Jacobian matrix onto the diagonal. v L = unknowns at the left and right boundaries. ˜ and add the transpose to obtain, We multiply (6) by v T H

I k ⊗ e 1T v

and v R =

˜ ˜J v + v T ˜J H˜ v t = B T L + B T R , v tT H

T I k ⊗ em v

are the

(7)

where

B T L = + v TL A − v L + g LT A + g L − ( v L − g L ) T A + ( v L − g L ), B T R = − v TR A + v R − g RT A − g R + ( v R − g R ) T A − ( v R − g R ).

(8)

Obtaining an energy estimate requires that the following holds

H J = J H,

(9)

˜ ˜J v t + v tT H˜ ˜J v = such that the LHS in (7) can be written as v H T

guaranteed if the norm H is diagonal.

d  v 2H J . dt

Hence, stability on curvilinear grids can only be

3. Computations in 1-D 3.1. Curvilinear grid generation Consider a smooth coordinate transformation x = x(ξ ), ξ ∈ [0, 1]. With the traditional (and extended) SBP operators, the grid points are equidistant in ξ , i.e., the traditional grid points are

ξi(T ) = (i − 1)h(T ) ,

i = 1, 2, . . . , m , h ( T ) =

1 m−1

where m denotes the number of grid points. On the stretched grid, the traditional grid points are (T )

xi

 = x ξi(T ) ,

i = 1, 2, . . . , m .

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

97

Table 2 log(l2 -errors) and convergence rates using 3 different 4th-order diagonal-norm SBP operators. The first (columns 2–3) shows the traditional, the second (columns 4–5) the extended, and the third (columns 6–7) the optimal. N

log l2

(T )

q( T )

log l2

(E)

q( E )

log l2

(O )

q( O )

51 101 201 401 801 1201

−1.09 −2.15 −3.12 −4.05 −4.96 −5.49

0.00 3.52 3.23 3.07 3.02 3.01

−1.11 −2.34 −3.56 −4.76 −5.94 −6.61

0.00 4.10 4.04 3.98 3.92 3.81

−1.25 −2.42 −3.61 −4.80 −6.00 −6.70

0.00 3.88 3.95 3.97 3.98 3.97

With the optimal SBP operators, the grid points are no longer equidistant in ξ ; the grid points are (O )

ξ1

= 0,

ξi( O ) = h( O )

i −1

dk ,

i = 2, 3, 4,

k =1

(O )

ξi

= h( O ) (d + i − 4),

ξm( O−)i = 1 − h( O )

i

dk ,

i = 5, 6, . . . , m − 4, i = 1, 2, 3,

k =1

(O )

ξm = 1, where h( O ) = 1/(2d + (m − 7)), d = d1 + d2 + d3 and hd1,2,3 are the grid-spacings to the right of the first, second, and third grid point, respectively. On the stretched grid, the grid points used by the optimal SBP operators are (O )

xi

(O ) , = x ξi

i = 1, 2, . . . , m .

3.2. First order systems We study the hyperbolic system (2) where





0 1 , A= 1 0

(10)

and where we specify the ingoing characteristics. Here



2A+ =



1 1 , 1 1



2A− =

−1 1



1 . −1

In this section we will make use of the following coordinate transformation,

 x = x(ξ ) = ξ

4 2

e −(ξ − 5 )

4 2

e −(1− 5 )

l ,

ξ ∈ [0, 1].

(11)

Here l is a non-negative integer, where the case l = 0 corresponds to a Cartesian grid. A larger value of l leads to a denser clustering of grid-points close to the boundaries. We approximate this system with the SBP-SAT method, given by (6). To test the accuracy of the newly developed diagonal-norm SBP operators we choose the analytic solution u = [sin(nπ (x − t )), − sin(nπ (x − t ))] T , where n is an integer number. In the present study we will set n = 10. The convergence rate is calculated as

 q = log10

  1/d m1 u − v (m1 ) h / log10 , m2 u − v (m2 ) h

(12)

where d is the dimension (d = 1 in the 1-D case), u is the analytic solution, and v (m1 ) the corresponding numerical solution with m1 unknowns. u − v (m1 ) h is the discrete l2 -norm of the error. In Tables 2–4 we compare the result using the traditional, extended and optimal diagonal-norm SBP operators, with no grid-stretching, i.e., l = 0 in (11). The solution is integrated in time to t = 0.4 using the 4th-order accurate Runge–Kutta method using a CFL of 0.5 (there is little difference between the time-step restrictions using the traditional and optimal SBP operators). In Fig. 2 we plot the l2 -errors (taken from Tables 3 and 4) versus the number of grid-points for the traditional and optimal 6th and 8th order accurate SBP operators in a log–log plot, to more clearly see the huge gain in efficiency going to optimal SBP operators.

98

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

Table 3 log(l2 -errors) and convergence rates using 3 different 6th order diagonal-norm SBP operators. The first (columns 2–3) shows the Traditional, the second (columns 4–5) the Extended, and the third (columns 6–7) the Optimal. N

log l2

(T )

q( T )

log l2

(E)

q( E )

log l2

(O )

q( O )

51 101 201 401 801 1201

−1.22 −2.52 −3.92 −5.35 −6.71 −7.47

0.00 4.31 4.65 4.75 4.54 4.29

−1.11 −2.71 −4.18 −5.68 −7.19 −8.07

0.00 5.32 4.88 4.98 5.00 5.00

−2.26 −4.04 −5.81 −7.53 −9.16 −10.08

0.00 5.91 5.88 5.71 5.42 5.24

Table 4 log(l2 -errors) and convergence rates comparing the traditional and optimal 8th-order diagonal-norm SBP operators. We also include the 10th order diagonal norm SBP operator with 11 boundary stencils. The first (columns 2–3) shows the Traditional 8th-order, the second (columns 4–5) the Optimal, 8th-order, and the third (columns 6–7) the 10th-order. N

log l2

(T )

q( T )

log l2

(O )

q( O )

log l2

(10th)

q(10th)

51 101 201 401 801 1201

−1.20 −2.68 −4.21 −5.76 −7.24 −8.04

0.00 4.92 5.10 5.13 4.91 4.58

−3.16 −5.41 −7.37 −9.04 −10.60 −11.46

0.00 7.47 6.49 5.55 5.20 4.89

−1.08 −2.73 −4.51 −6.31 −8.14 −9.21

0.00 5.49 5.92 5.97 6.08 6.09

Fig. 2. Comparing the traditional and optimal 6th- and 8th-order accurate SBP operators in a log–log plot, results taken from Tables 3 and 4.

To validate the importance of using diagonal-norm SBP operators for this problem we present in Fig. 3 the eigenvalues to (6) where the system is defined by (10), on a grid defined by (11) with l = 2 and m = 51. We compare the spectra using the 8th order optimal diagonal-norm SBP operator and the corresponding block-norm SBP operator. The largest real part of the spectrum is referred to as λmax . A non-positive value is required for stability. 3.3. Second order systems We study the second order wave equation on piecewise continuous media, (l)



(l) 

a(l) utt = b(l) u x a

(r ) (r )



utt = b (l)

x

ux

(r )

b(l) u x = b(r ) u x ,

−1  x  0, t  0,

,

 (r ) (r )

x

,

0  x  1, t  0, u (l) = u (r ) ,

x = 0, t  0,

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

99

Fig. 3. The spectrum to (6), where the system is defined by (10), comparing the 8th order optimal diagonal-norm SBP operator and the 8th order block-norm SBP operator. m = 51. λmax represents the largest real part of the spectrum.

(l)

u x = 0,

x = −1, t  0,

(r )

u x = 0,

x = 1 , t  0,

u

(l)

=f ,

u

(r )

(r )

(l)

(l)

=f

u t = 0, (r )

ut

,

−1  x  0, t = 0,

= 0,

0  x  1 , t = 0,

(13)

where a(l) = a(r ) , b(l) = b(r ) . Here u (l,r ) denote the solutions corresponding to the left and right domains respectively, and f (l,r ) the corresponding initial data. (A detailed analysis of this particular problem is described in [29].) We introduce the following coordinate transformations,

 x = x(ξ ) = ξ

l ,

4 2

 x = x(ξ ) = ξ

4 2

e −(ξ + 5 )

e −(−1+ 5 ) 4 2

e −(ξ − 5 )

ξ ∈ [−1, 0],

l ,

4 2

e −(1− 5 )

ξ ∈ [0, 1],

(14)

in the left and right domains, respectively. Here l is a non-negative integer, where the case l = 0 corresponds to the Cartesian case. A larger value of l leads to a denser clustering of grid-points close to the interface and the outer boundaries. The problem (13) transforms to





(l) (l) a˜ (l) utt = b˜ (l) u ξ ξ ,



(r ) 

(r ) a˜ (r ) utt = b˜ (r ) u ξ

−1  ξ  0, t  0, ,

0  ξ  1, t  0,

ξ (l) (r ) b˜ (l) u ξ = b˜ (r ) u ξ , u (l) = u (r ) , (l) u ξ = 0, ξ = −1, t  0, (r ) u ξ = 0, ξ = 1, t  0, (l) (l) (l)

=f ,

u u

(r )

=f

(r )

,

u t = 0, (r )

ut

= 0,

ξ = 0, t  0,

−1  ξ  0, t = 0, 0  ξ  1 , t = 0,

(15)

(l,r ) (l,r ) in curvilinear coordinates, where a˜ (l,r ) = a(l,r ) xξ and b˜ (l,r ) = b(l,r ) /xξ .

The semi-discrete approximation of the continuity conditions u (l) = u (r ) and b˜ (l) u ξ = b˜ (r ) u ξ (l)

T (l) I 1 ≡ em v − e 1T v (r ) = 0,

I 2 ≡ B˜ (l) em D 1 v (l) − B˜ (r ) e 1 D 1 v (r ) = 0,

(r )

can be written,

(16)

where v (l,r ) are the solution vectors corresponding to the left and right domains respectively, and B˜ (l,r ) are the projections of b˜ l,r (ξ ) onto the diagonals. The left and right domains are discretized using m grid points.

100

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

Fig. 4. The spectrum to (17) on a grid defined by (14) with l = 2 and m = 51. Comparing the optimal diagonal 8th order SBP operator and the 8th order block-norm SBP operator. Imaginary eigenvalues lead to instability.

The SBP-SAT method for (15) is given by,

˜ (l) v tt = D 1 B˜ (l) D 1 v (l) + τ H −1 em ( I 1 ) + β H −1 B˜ (l) (em D 1 )T ( I 1 ) + γ H −1 em ( I 2 ) − H −1 B˜ (l) e 1 D 1 v (l) , A (l)

˜ (r ) v tt = D 1 B˜ (r ) D 1 v (r ) − τ H −1 e 1 ( I 1 ) − β H −1 B˜ (l) (e 1 D 1 )T ( I 1 ) − γ H −1 e 1 ( I 2 ) + H −1 B˜ (r ) em D 1 v (r ) . A (2 )

(17)

The following lemma was first stated in [29]: ˜ (l) ˜ (r )

b Lemma 3.1. The scheme (17) is strictly stable if D 1 is a diagonal-norm SBP operator and γ = − 12 , β = 12 , and τ  − b 4/+5h

hold.

The proof can be found in [29]. To show the importance of using diagonal-norm SBP operators for this problem we present in Fig. 4 the eigenvalues to (17), on a grid defined by (14) with l = 2 and m = 51. We compare the spectra using the optimal 8th order diagonal-norm SBP operator and the corresponding block-norm SBP operator. Stability cannot be guaranteed on curvilinear grids using block-norm SBP operators. (Stability requires that all eigenvalues are real and non-positive.) The accuracy of the discontinuous media treatment with the SBP-SAT method, given by (17) will be verified against an analytic solution. A narrow Gaussian pulse is propagated across the discontinuous media interface to verify the accuracy of the transmission and reflection properties at the interface. The initial data is given by



u (x, 0) = exp −

(x0 − x)2 r02



,

ut (x, 0) = 0,

where x0 determines the position and r0 determines the width of the Gaussian profile. In the present study we set x0 = − 14 and r0 =

1 . 30

Hence, the initial Gaussian profile is centred at x = −0.25 in the left sub-domain where the media parameters

(l) (l) are given  by a , b . The speed of sound and the acoustic impedance in the left and right subdomains are given by

c (l,r ) =

a(l,r ) , b(l,r )

√

σ (l,r ) = a(l,r ) b(l,r ) . For t > 0 the initial Gaussian profile splits into one left-going and one right-going Gaussian

wave (see Fig. 5) propagating with speed c (l) , i.e., the analytic solution is given by

u (x, t ) =

1 2



exp −

(x0 − (x + c (l)t ))2



r02

+

1 2



exp −

(x0 − (x − c (l)t ))2 r02



, −x

assuming that the waves have not yet reached the left boundary or the interface. At t ∗ = (l)0 (in the present study t ∗ = c 0.25), the right-going wave hits the media interface at x = 0 and splits into a reflected wave and a transmitted wave, see Fig. 5. The transmission (T ) and reflection (R) coefficients are given by

T=

2σ (l)

σ (l)

+ σ (r )

,

R=

σ (l) − σ (r ) . σ (l) + σ (r )

Hence, at t > t ∗ the analytic solution consists of a left-going wave (with speed c (l) and amplitude condition, a reflected wave with amplitude transmitted right-going wave with amplitude

R 2 T 2

1 ) 2

from the initial

moving (with speed c (l) ) to the left from the interface at x propagating with speed c (r ) . The analytic solution is given by

= 0, and a

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

101

Fig. 5. The numerical solution at t = 0.15, t = 0.25, t = 0.30 and t = 0.50, solving (17) on a grid defined by (14) with l = 0 and m = 101, using the optimal 8th order diagonal-norm SBP operator, corresponding to the results using m = 101 in Table 6.

u

(l)

=

1 2

    (x0 − (x + c (l)t ))2 (x + c (l) (t − t ∗ ))2 R + exp − x ∈ [−1, 0], exp − 2 2 

r0

2

 (c (l) (x − c (r ) (t − t ∗ )))2 T , u (r ) = exp − 2 (c (r ) r0 )2

r0

x ∈ [0, 1],

(18)

assuming that left-going and right-going waves have not yet reached the outer boundaries. In the numerical simulations we choose a(l) = 1, b(l) = 1, a(r ) = 2, b(r ) = 8. At the outer boundaries we impose homogeneous Neumann boundary conditions. The numerical approximations are integrated to t = 0.5 using the 4th order Runge–Kutta method, with a time step dt = 0.05h. (There is little difference between the time-step restrictions using the traditional and optimal SBP operators.) In Tables 5 and 6 we compare the result using the traditional and optimal diagonalnorm SBP operators. We compare the convergence on both an equidistant grid, by setting l = 0 in (14), and a stretched grid, by setting l = 2 in (14). In Fig. 6 we plot the numerical and analytical solution using the two different 8th order SBP operators, corresponding to the results using m = 101 in Table 6. 4. The Euler vortex problem To test the accuracy properties of the new optimal diagonal-norm SBP operators in a more challenging 2-D application, we consider the nonlinear Euler equations (see for example [40]) in curvilinear coordinates formulated in conservative form

102

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

Fig. 6. The numerical (solid line) and analytical (dotted line) solution at t = 0.5, solving (17) on a grid defined by (14) with l = 0 and m = 101. Comparing the traditional (left) and optimal (right) 8th order diagonal-norm SBP operators.

Table 5 log(l2 -errors) and convergence rates, discontinuous interface. Comparing the Traditional (columns 2–3) and Optimal (columns 4–5) 6th order diagonal-norm SBP operators on non-stretched grids, l = 0 in (14). The last four columns show the traditional and optimal 6th order diagonal-norm on a stretched grid, l = 2 in (14). (T )

(O )

(T )

(O )

N

log l2 (l=0)

q

log l2 (l=0)

q

log l2 (l=2)

q

log l2 (l=2)

q

51 101 201 401 801

−1.42 −2.22 −3.09 −4.36 −5.75

0.00 2.69 2.92 4.23 4.63

−1.45 −2.58 −4.28 −6.06 −7.85

0.00 3.80 5.69 5.94 5.95

−1.27 −2.06 −3.55 −5.31 −7.05

0.00 2.65 4.99 5.87 5.79

−1.33 −2.05 −3.54 −5.31 −7.05

0.00 2.43 4.99 5.87 5.79

Table 6 log(l2 -errors) and convergence rates, discontinuous interface. Comparing the Traditional (columns 2–3) and Optimal (columns 4–5) 8th order diagonal-norm SBP operators on non-stretched grids, l = 0 in (14). The last four columns show the traditional and optimal 8th order diagonal-norm on a stretched grid, l = 2 in (14). (T )

(O )

(T )

(O )

N

log l2 (l=0)

q

log l2 (l=0)

q

log l2 (l=2)

q

log l2 (l=2)

q

51 101 201 401 801

−1.44 −2.20 −3.76 −4.99 −6.48

0.00 2.58 5.22 4.09 4.97

−1.65 −3.12 −5.38 −7.75 −9.81

0.00 4.94 7.58 7.89 6.85

−1.45 −2.38 −4.38 −6.70 −8.71

0.00 3.15 6.68 7.75 6.67

−1.40 −2.36 −4.35 −6.70 −9.09

0.00 3.24 6.68 7.82 7.97

˜ η = 0, J ut + F˜ ξ + G

(19)

where

⎡

ρ ⎤ ⎢ ρu ⎥ u=⎣ , ρv ⎦

⎡

ρU

⎤

⎢ ρ uU + y η p ⎥ , ρ vU − xη p ⎦ (e + p )U

F˜ = ⎣

e

⎡

ρV

⎤

⎢ ρ uV − y ξ p ⎥ . ρ v V + xξ p ⎦ (e + p ) V

˜ =⎣ G

Here

J = xξ y η − xη y ξ , U = y η u − xη v ,

V = − y ξ + xξ v .

ρ is the density (scaled to free stream density ρ∞ ); u, v are the velocity components in the x, y directions respectively ρ 2 ) (scaled to free stream speed of sound c ∞ ); e = e i + 2 (u 2 + v 2 ) the total energy (nondimensionalized with respect to ρ∞ c ∞ and e i is the internal energy. Pressure is obtained from the equation of state for a perfect gas p = (γ − 1)e i .

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

103

In the application below we consider the propagation of a vortex that satisfies the two-dimensional Euler equations, under the assumption of isentropy. The analytic vortex solution is steady in the frame of reference moving with the freestream. The scaled vortex has the velocity field



vΘ =



r 1 − r2 , exp 2π 2

(20)

where  is the nondimensional circulation, (r , Θ) are the polar coordinates. The analytic solution in a fixed frame of reference (x, y , t ) becomes



u=1−



y f (x, y , t ) , exp 2π 2



v=





 ((x − x0 ) − t ) f (x, y , t ) , exp 2π 2  γ −1 1

  2 (γ − 1) M 2 ρ = 1− exp f (x, y , t ) 2 8π

,

p=

ργ , γ M2

(21)

where f (x, y , t ) = 1 − (((x − x0 ) − t )2 + y 2 ); M is the Mach number; γ = c p /c v = 1.4 and x0 is the initial position of the vortex (in the x-direction). The vortex is introduced into the computational domain by using the analytic solution as boundary data and initial data. The main focus is to verify the accuracy for a non-smooth multi-block interface coupling. Hence, we will solve the Euler-vortex problems on a two-block domain (see Figs. 8 and 9). Details concerning how to discretize this problem with the SBP-SAT method can be found in [47,33,27]. 4.1. Curvilinear grid generation in 2-D Consider a smooth one-to-one mapping from the unit square to a curvilinear domain Ω ,



x = x(ξ, η) y = y (ξ, η)

(ξ, η) ∈ [0, 1] × [0, 1], (ξ, η) ∈ [0, 1] × [0, 1].

(22)

The traditional grid points in the unit square are





 ξi(T ) , η(jT ) = (i − 1)hξ(T ) , ( j − 1)h(ηT ) ,

i = 1, . . . , m ξ , j = 1, . . . , m η ,

(T )

where hξ,η = 1/(mξ,η − 1) and mξ,η are the numbers of grid points in the ξ - and grid points in the curvilinear domain Ω are



(T )

(T ) 

xi j , y i j





 = x ξi(T ) , η(jT ) , y ξi(T ) , η(jT ) ,

η-directions, respectively. The traditional

i = 1, . . . , m ξ , j = 1, . . . , m η .

The grid points in the unit square associated with the optimal SBP operators are



 ξi( O ) , η(j O ) ,

i = 1, . . . , m ξ , j = 1, . . . , m η ,

where

ξ1( O ) = 0, ξi( O ) = h(ξO )

η1( O ) = 0, i −1 k =1

dk ,

i = 2, 3, 4,

ξi( O ) = h(ξO ) (d + i − 4), (O )

(O )

ξmξ −i = 1 − hξ

i  k =1

i = 5, . . . , m ξ − 4,

dk ,

i = 1, 2, 3,

(O )

η(j O ) = h(ηO )

j −1 k =1

dk ,

j = 2, 3, 4,

η(j O ) = h(ηO ) (d + j − 4), (O ) (O ) ηm = 1 − hη η− j

j  k =1

dk ,

j = 5, . . . , m η − 4, j = 1, 2, 3,

(O ) ηm η = 1,

ξmξ = 1, and

1

(O )

hξ,η =

2d + (mξ,η − 7)

.

Thus, the non-equal spacings between the grid points appear in each dimension separately. This makes the extension to higher dimensions straightforward. The grid points used by the optimal SBP operators in the curvilinear domain Ω are



(O )

(O )

xi j , y i j





 = x ξi( O ) , η(j O ) , y ξi( O ) , η(j O ) ,

i = 1, . . . , m ξ , j = 1, . . . , m η .

104

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

Fig. 7. Comparing the traditional and optimal 6th and 8th order accurate SBP operators in a log–log plot, results taken from Tables 7 and 8.

4.2. Computations In the first two tests we verify the convergence properties using non-stretched grids (see Fig. 8) and stretched grids (see Fig. 9). The initial position of the centre of the vortex is (x, y ) = (0, −0.5), so that the grid-interface—where there is a discontinuity in the metrics—cuts the vortex (see Fig. 10). The vortex is then propagated (to the right) along the grid interface a distance of 4 units (corresponding to t = 4) where we measure the l2 -errors. The vortex is integrated using a 6th-order accurate Runge–Kutta method [2] using a CFL of 0.4. We set the Mach number M = 0.5. In Tables 7 and 8 we present the results using the traditional diagonal-norm and optimal diagonal-norm SBP operators where we apply no stretching, corresponding to l = 0 in (14) in both the x- and y-direction, see Fig. 8. We performed a detailed study (although not presented here) of what would be the optimal (in terms of accuracy) grid-stretching for the traditional 8th order diagonal SBP operator for this 2-block grid. This study yielded a stretching corresponding to l = 2 in (14) in both the x- and y-direction. In Tables 10 and 11 we present the results using the traditional diagonal-norm and optimal diagonal-norm SBP operators where we apply the optimal grid-stretching corresponding to l = 2 in (14) in both the x- and y-direction, see Fig. 8. In Table 9 we list the CFL constraints (relative the 2nd-order accurate traditional SBP scheme on a non-stretched grid) for the traditional and optimal SBP schemes for the Euler vortex problem, for both the non-stretched grid, l = 0 in (14) and the stretched grid, l = 2 in (14). In Fig. 7 we plot the l2 -errors (taken from Tables 7 and 8) versus the number of grid-points (in each dimension) for the traditional and optimal 6th- and 8th-order accurate SBP operators in a log–log plot, to more clearly see the huge gain in efficiency going to optimal SBP operators. Remark. By comparing Tables 8 and 10 we see that the optimal SBP operators yield accurate solutions on coarse grids, without grid-stretching. By comparing Tables 8 and 11 it is clear that stretching does not improve accuracy for the optimal SBP. Grid-stretching (here corresponding to l = 2 in (14)) does however improve accuracy for the traditional SBP. We see from Table 9 that we get a relative CFL of 0.37 for the optimal 8th-order SBP with no stretching and 0.21 for the traditional 8th-order SBP with stretching. Hence we gain almost a factor 2 in favour of the optimal SBP (assuming we use l = 0 for the optimal and l = 2 for the traditional) concerning the time-step restriction. A fair comparison of efficiency is given by analysing Tables 8 and 10. We need approximately twice as fine grid using the traditional 8th-order SBP (on a stretched grid) to obtain the same error as with the 8th-order optimal (non-stretched grid). Here we gain 2d , where d is the number of spatial dimensions. We then roughly get an additional factor of 4 taking the time-step into consideration, i.e., the gain (in floating point operations), in favour of the optimal 8th-order SBP would be roughly 4 × 23 = 32, for a 3-D problem. To test the stability properties and also the superior accuracy of the optimal diagonal-norm SBP operators (compared to the traditional diagonal-norm SBP operators) we perform long-time simulations to t = 34. The west and east boundaries are coupled using SAT, and it takes 12 units in time for the vortex to cross the domain (from left to right). At t = 34 the vortex has crossed the west-east interface 3 times.

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

105

Fig. 8. The computational domain, consisting of 2 curvilinear grids, with discontinuous metrics along the interface. Here with 21 × 11 grid-points in each domain and no stretching of the grids, l = 0 in (14).

Fig. 9. The computational domain, consisting of 2 curvilinear grids, with discontinuous metrics along the interface. Here with 21 × 11 grid-points in each domain and stretched grids, l = 2 in (14).

Table 7 log(l2 -errors) and convergence rates using the traditional diagonal-norm operators with no grid stretching l = 0. N x and N y denote the number of grid points in the x and y directions, respectively, in each subdomain. Nx × N y 41 × 21 81 × 41 161 × 81 321 × 161

(2nd)

(4th)

(6th)

(8th)

log l2

q(2nd)

log l2

q(4th)

log l2

q(6th)

log l2

q(8th)

−2.41 −3.06 −3.66 −4.27

0.00 2.17 2.00 2.00

−2.70 −3.69 −4.67 −5.60

0.00 3.29 3.25 3.09

−2.23 −3.38 −4.75 −6.29

0.00 3.81 4.58 5.09

−2.40 −3.49 −4.88 −6.33

0.00 3.60 4.62 4.82

In Fig. 11 we present the solutions for two different schemes: traditional 8th order diagonal-norm SBP and optimal 8th order diagonal-norm SBP, using 31 × 61 grid-points in each subdomain. The optimal 8th order diagonal-norm SBP operator yields much more accurate results. Both SBP operators are closed with 4th order accurate boundary stencils, and it is evident that the accuracy of the boundary stencils dictates the accuracy properties when block-interfaces (or boundaries) are present.

106

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

Fig. 10. The solution (pressure) at t = 0 (a), and later at t = 6 (b) when the vortex reaches the outermost right boundary, and re-enters through the left boundary.

Fig. 11. The solution (pressure) at t = 34 using m = 31 × 61 grid-points (in each sub-block), comparing the solutions with the traditional and optimal 8th order SBP.

Table 8 log(l2 -errors) and convergence rates using the optimal diagonal-norm operators with no grid stretching l = 0. N x and N y denote the number of grid points in the x and y directions, respectively, in each subdomain. Nx × N y 41 × 21 81 × 41 161 × 81 321 × 161

(2nd)

(4th)

(6th)

(8th)

log l2

q(2nd)

log l2

q(4th)

log l2

q(6th)

log l2

q(8th)

−2.51 −3.12 −3.71 −4.31

0.00 2.02 1.98 1.99

−3.45 −4.63 −5.72 −6.84

0.00 3.91 3.63 3.71

−3.49 −5.09 −6.57 −7.98

0.00 5.33 4.91 4.66

−3.75 −5.82 −7.91 −9.84

0.00 6.89 6.93 6.43

5. Conclusions and future work The main focus has been to construct new diagonal-norm SBP operators with optimal accuracy properties, to combine truly high-order accuracy and stability on curvilinear multi-block grids. This is achieved by introducing more boundary stencils and (or) allowing the grid-points to be non-equidistant close to the boundaries. The boundary and interface conditions

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

107

Table 9 Maximal time step with no grid stretching (l = 0) and with l = 2 in (14), relative the time-step for the 2nd order scheme on a non-stretched grid. Determined experimentally on an 81 × 41 grid. Operator

2ndl=0

4thl=0

6thl=0

8thl=0

2ndl=2

4thl=2

6thl=2

8thl=2

Traditional Optimal

1 0.93

0.60 0.58

0.55 0.49

0.49 0.37

0.45 0.28

0.27 0.14

0.22 0.12

0.21 0.09

Table 10 log(l2 -errors) and convergence rates using the traditional diagonal-norm operators with grid stretching l = 2. N x and N y denote the number of grid points in the x and y directions, respectively, in each subdomain. (2nd)

Nx × N y

(4th)

(6th)

(8th)

log l2

q(2nd)

log l2

q(4th)

log l2

q(6th)

log l2

q(8th)

−2.24 −2.79 −3.38 −3.98

0.00 1.80 1.98 2.00

−2.91 −4.00 −5.14 −6.27

0.00 3.60 3.81 3.76

−2.84 −4.47 −6.15 −7.73

0.00 5.42 5.58 5.22

−2.57 −4.17 −5.93 −7.33

0.00 5.29 5.85 4.65

41 × 21 81 × 41 161 × 81 321 × 161

Table 11 log(l2 -errors) and convergence rates using the optimal diagonal-norm operators with grid stretching l = 2. N x and N y denote the number of grid points in the x and y directions, respectively, in each subdomain. (2nd)

Nx × N y

(4th)

(6th)

(8th)

log l2

q(2nd)

log l2

q(4th)

log l2

q(6th)

log l2

q(8th)

−2.24 −2.78 −3.38 −3.98

0.00 1.81 1.99 2.00

−2.93 −4.02 −5.18 −6.38

0.00 3.59 3.87 3.97

−3.27 −4.88 −6.64 −8.43

0.00 5.34 5.85 5.96

−3.46 −5.62 −7.97 −10.36

0.00 7.17 7.80 7.97

41 × 21 81 × 41 161 × 81 321 × 161

are treated with the SAT technique. We study both first and second order hyperbolic problems on piecewise curvilinear multi-block grids. Numerical computations in 1-D and 2-D corroborate the stability- and accuracy-properties and show that the optimal diagonal-norm SBP operators are superior to the corresponding traditional diagonal-norm SBP operators. For the second order wave equation we obtain super-convergence, and also for the 2-D Euler vortex problem the new optimal SBP operators yield higher than expected order of convergence. It is likely that further improvements in accuracy could have been achieved for the optimal 6th and 8th-order cases by allowing for more than 3 “free” locations of the boundary points in the minimization. However, due to a limitation in Maple we where forced to limit the number to 3. With other minimization routines this limitation might be lifted. The 8th-order case could potentially allow for 7 free boundary points. In a coming study we will investigate the efficiency of using the new diagonal-norm SBP operators in a more complex 3-D setting, with focus on second order hyperbolic systems and hyperbolic–parabolic systems. Appendix A. Extended diagonal-norm SBP operators The boundary closures (the coefficients) are presented for the Q and H matrices. Notice that the coefficients in the norm should be multiplied with the grid-step h to have the correct normalisation. A.1. 2nd-order case The upper part of the norm H

H 1,1 =

5 12

H 2,2 =

,

7 6

H 3,3 =

,

11 12

.

The upper part of Q ,

Q 1,2 =

7 12

Q 1,3 = −

,

1 12

Q 2,3 =

,

7 12

.

A.2. 4th-order case The upper part of the norm H ,

H 1,1 = H 2,2 =

511 1600 3971 2880

,

H 3,3 =

,

H 4,4 =

929 1440 587 480

,

,

H 5,5 = H 6,6 =

2659

,

2880 14 551 14 400

.

108

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

The upper part of Q ,

Q 1,2 =

436 061

Q 1,3 = − Q 1,4 = − Q 1,5 = Q 1,6 =

Q 2,3 =

,

680 400 312 391 5 443 200 100 109

907 200 54 107

5 443 200 263 15 552

Q 2,4 =

,

238 831

1 088 640 125

Q 2,5 = −

,

Q 2,6 = −

,

Q 3,4 =

,

,

544 320 290 357

Q 3,5 = −

,

Q 3,6 =

,

5184 220 357

5 443 200 3467 7776

Q 4,5 =

,

74 827

1 088 640 3859

907 200 342 799

,

,

544 320 149 959

Q 4,6 = −

,

Q 5,6 =

5 443 200 428 789 680 400

,

,

.

A.3. 6th-order case The upper part of the norm H ,

H 1,1 = H 2,2 =

7 497 391 25 401 600 1 106 227 725 760

H 3,3 =

, H 4,4 =

,

H 5,5 =

105 317 403 200 260 179 145 152 303 631 725 760

,

H 6,6 =

,

H 7,7 =

,

H 8,8 =

513 973 403 200 671 171

, ,

725 760 25 631 999 25 401 600

.

Left boundary closure of Q ,

Q 1,2 =

14 085 349 21 168 000 205 189

Q 1,3 = − Q 1,4 = − Q 1,5 = − Q 1,6 =

80

Q 1,8 = − Q 2,3 = Q 2,4 =

38 102 400 89 573

Q 2,7 =

,

Q 2,8 =

,

Q 3,4 =

,

39 279 943 1 524 096 000 13

,

2000 56 626 691

304 819 200 97 579 603 152 409 600

1 214 239

,

16 934 400 246 914 711

Q 2,6 = −

,

8 064 000 7 583 609

4 233 600 9

Q 1,7 = −

Q 2,5 =

,

Q 4,5 =

, ,

,

762 048 000 23 380 423

304 819 200 3 200 22 302 157

Q 5,6 =

101 606 400 6 048 877

Q 5,7 =

60 963 840 9 737 351

169 344 000 508 309

,

304 819 200 2 075 677

152 409 600 29 205 107

152 409 600 15 263 287 304 819 200 81 127

Q 6,7 =

,

,

508 032 000 49 594 423

Q 6,8 = −

, Q 7,8 =

762 048 000 30 841 499 43 545 600

,

,

,

3 386 880 321 139 631

,

,

Q 5,8 = −

,

25 401 600 780 937

304 819 200

,

Q 4,8 = −

,

,

6 350 400 42 123 859

Q 4,7 = −

,

Q 3,7 = − Q 3,8 =

,

30 481 920 3 467 971

,

Q 3,5 = − Q 3,6 =

Q 4,6 =

8 120 941

,

.

Appendix B. Optimal diagonal-norm SBP operators, with three non-equidistant grid-points Let x ∈ [0, 1]. Introduce d = d1 + d2 + d3 , where hd1 is the distance between the first and second grid-point, hd2 is the distance between the second and third grid-point, and hd3 is the distance between the third and fourth grid-point. Here h = 1/(2d + (m − 7)) and m is the number of grid-points. The grid-points are located at









x = 0, d1 h, (d1 + d2 )h, dh, (d + 1)h, (d + 2)h, . . . , d + (m − 7) h, 1 − dh, 1 − (d1 + d2 )h, 1 − d1 h, 1 .

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

B.1. 2nd-order case The first three grid-stretchings,

d1 = 0.78866488858096586513,

d2 = 0.95915098594220826013,

d 3 = 1.

The upper part of the norm H ,

H 1,1 = 0.33743097329453577701, H 2,2 = 0.97759682018833491296,

H 3,3 = 0.93278808104030343530.

The upper part of the Q ,

Q 1,2 = 0.55932483188770411252, Q 1,3 = −0.05932483188770411252,

Q 2,3 = 0.55932483188770411252.

B.2. 4th-order case The first three grid-stretchings,

d1 = 0.72181367003646814327, d2 = 1.3409118421582217252,

d3 = 1.2898797485951900258.

The upper part of the norm H ,

H 3,3 = 1.434458792494126, H 4,4 = 1.0917323021736130836, H 5,5 = 0.9883816115129601975.

H 1,1 = 0.21427296612044126417, H 2,2 = 1.123759588488739348, The upper part of Q ,

Q 1,2 = 0.66884898686930380508, Q 1,3 = −0.25171531878753856238, Q 1,4 = 0.10997619816825822803, Q 1,5 = −0.027109866250023470592, Q 2,3 = 0.92214436948640491071,

Q 2,4 = −0.32412368653542520402, Q 2,5 = 0.070828303918324098284, Q 3,4 = 0.8180378089216779335, Q 3,5 = −0.14760875822281158529, Q 4,5 = 0.68722365388784429092.

B.3. 6th-order case The first three grid-stretchings,

d1 = 0.51670081689316731234, d2 = 0.98190527037374634269,

d3 = 1.0868393364992957832.

The upper part of the norm H ,

H 1,1 = 0.15109714532036117328, H 2,2 = 0.80967585357107013003, H 3,3 = 1.0911427148079254850,

H 4,4 = 1.0435269041571577756, H 5,5 = 0.98680905919946100728, H 6,6 = 1.0037581831426163456, H 7,7 = 0.99943556356761752125.

The upper part of Q ,

Q 1,2 = 0.66670790901888837033, Q 1,3 = −0.23418791580399147484, Q 1,4 = 0.084251264588860596867, Q 1,5 = −0.015923290838179674350, Q 1,6 = −0.0015653772860347171721, Q 1,7 = 0.00071741032045689717567, Q 2,3 = 0.89405599296515541581, Q 2,4 = −0.28597427787314667440, Q 2,5 = 0.057056178538117177397, Q 2,6 = 0.0041320613074890940489, Q 2,7 = −0.0025620459187266476645,

Q 3,4 = 0.82961715259707113283, Q 3,5 = −0.18233747042994439227, Q 3,6 = 0.0083784382166533084621, Q 3,7 = 0.0042099567773838744673, Q 4,5 = 0.75419218459746682761, Q 4,6 = −0.14034899831339963049, Q 4,7 = 0.014050953028717865444, Q 5,6 = 0.74751473989919011204, Q 5,7 = −0.15119380469839682133, Q 6,7 = 0.75144419715723149872.

109

110

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

B.4. 8th-order case The first three grid-stretchings

d1 = 0.41669687672575697416, d2 = 0.78703773886730090312,

d3 = 0.92685925671601406028.

The upper part of the norm H ,

H 1,1 = 0.12163222110707502878, H 2,2 = 0.65235832636546639982, H 3,3 = 0.87730414198101010954, H 4,4 = 0.97388951771079542799,

H 5,5 = 1.0072514376844677230, H 6,6 = 0.99768726657776478834, H 7,7 = 1.0005302998791085514, H 8,8 = 0.99994066100338390832.

The upper part of Q ,

Q 1,2 = 0.66670790901888837033, Q 1,3 = −0.21994030190635039046, Q 1,4 = 0.061752567584332553851, Q 1,5 = −0.0032312350944133128873, Q 1,6 = −0.0033934980320003350186, Q 1,7 = 0.000015157027970563223705, Q 1,8 = 0.00032350133072942893419, Q 2,3 = 0.86688767821045233147, Q 2,4 = −0.24298087640343350527, Q 2,5 = 0.039549469619698650847, Q 2,6 = 0.0020763528371484737510, Q 2,7 = −0.00065045489396961912976, Q 2,8 = −0.00040836028016484152445,

Q 3,4 = 0.82065092584472835146, Q 3,5 = −0.21014872891771196683, Q 3,6 = 0.038572177503610408523, Q 3,7 = −0.0016637807199883547459, Q 3,8 = −0.00046321740653650899692, Q 4,5 = 0.81102837324727866266, Q 4,6 = −0.20423896795865859484, Q 4,7 = 0.034947121761434371005, Q 4,8 = 0.0023139100244270367378, Q 5,6 = 0.80054065093594950025, Q 5,7 = −0.19699167992690472530, Q 5,8 = 0.037220336417235830241, Q 6,7 = 0.79903079313046586260, Q 6,8 = −0.19999788736822592697, Q 7,8 = 0.80016334685519857774.

References [1] S. Abarbanel, A. Ditkowski, Asymptotically stable fourth-order accurate schemes for the diffusion equation on complex shapes, J. Comput. Phys. 133 (1997) 279–288. [2] E.A. Alshinaa, E.M. Zaksb, N.N. Kalitkin, Optimal first- to sixth-order accurate Runge–Kutta schemes, Comput. Math. Math. Phys. 48 (2008) 395–405. [3] A. Bayliss, K.E. Jordan, B.J. Lemesurier, E. Turkel, A fourth order accurate finite difference scheme for the computation of elastic waves, Bull. Seismol. Soc. Am. 76 (4) (1986) 1115–1132. [4] M.H. Carpenter, D. Gottlieb, S. Abarbanel, The stability of numerical boundary treatments for compact high-order finite difference schemes, J. Comput. Phys. 108 (2) (1993). [5] M.H. Carpenter, D. Gottlieb, S. Abarbanel, Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes, J. Comput. Phys. 111 (2) (1994). [6] M.H. Carpenter, J. Nordström, D. Gottlieb, A stable and conservative interface treatment of arbitrary spatial accuracy, J. Comput. Phys. 148 (1999). [7] Mark H. Carpenter, Travis C. Fisher, High-Order Entropy Stable Formulations for Computational Fluid Dynamics, American Institute of Aeronautics and Astronautics, 2013, 2013/06/27. [8] Mark H. Carpenter, David Gottlieb, Spectral methods on arbitrary grids, J. Comput. Phys. 129 (1996) 74–86. [9] Peter Diener, Ernst Nils Dorband, Erik Schnetter, Manuel 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) (2007) 109–145. [10] Travis H. Fisher, Mark H. Carpenter, High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains, J. Comput. Phys. 252 (2013) 518–557. [11] D. Funaro, D. Gottlieb, A new method of imposing boundary conditions in pseudospectral approximations of hyperbolic equations, Math. Comput. 51 (184) (1988) 599–613. [12] D. Funaro, D. Gottlieb, Convergence results for pseudospectral approximations of hyperbolic systems by a penalty-type boundary treatment, Math. Comput. 57 (196) (1991) 585–596. [13] B. Gustafsson, The convergence rate for difference approximations to general mixed initial boundary value problems, SIAM J. Numer. Anal. 18 (2) (Apr. 1981) 179–190. [14] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time Dependent Problems and Difference Methods, John Wiley & Sons, Inc., 1995. [15] B. Gustafsson, H.O. Kreiss, A. Sundström, Stability theory of difference approximations for mixed initial boundary value problems, Math. Comput. 26 (119) (1972). [16] B. Gustafsson, P. Olsson, Fourth-order difference methods for hyperbolic IBVPs, J. Comput. Phys. 117 (1) (1995). [17] J.S. Hesthaven, A stable penalty method for the compressible Navier–Stokes equations: II Multi-dimensional domain decomposition schemes, SIAM J. Sci. Comput. (1998). [18] J. Hicken, D. Zingg, Superconvergent functional estimates from summation-by-parts finite-difference discretizations, SIAM J. Sci. Comput. 33 (2) (2011) 893–922. [19] J.E. Hicken, D.W. Zingg, Parallel Newton–Krylov solver for the Euler equations discretized using simultaneous approximation terms, AIAA J. 46 (11) (2008) 2773–2786, 2013/02/09.

K. Mattsson et al. / Journal of Computational Physics 264 (2014) 91–111

111

[20] Jason E. Hicken, David W. Zingg, Aerodynamic optimization algorithm with integrated geometry parameterization and mesh movement, AIAA J. 48 (2) (2010) 400–413, 2013/02/09. [21] J.E. Hicken, Output error estimation for summation-by-parts finite-difference schemes, J. Comput. Phys. 231 (9) (2012) 3828–3848. [22] H.-O. Kreiss, G. Scherer, Finite element and finite difference methods for hyperbolic partial differential equations, in: Mathematical Aspects of Finite Elements in Partial Differential Equations, Academic Press, Inc., 1974. [23] H.-O. Kreiss, G. Scherer, On the existence of energy estimates for difference approximations for hyperbolic systems, Technical report, Dept. of Scientific Computing, Uppsala University, 1977. [24] Heinz-Otto Kreiss, Joseph Oliger, Comparison of accurate methods for the integration of hyperbolic equations, Tellus XXIV (3) (1972). [25] L. Lehner, O. Reula, M. Tiglio, Multi-block simulations in general relativity: high-order discretizations, numerical stability and applications, Class. Quantum Gravity 22 (2005) 5283–5321. [26] S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42. [27] K. Mattsson, M. Almquist, A solution to the stability issues with block norm summation by parts operators, J. Comput. Phys. 253 (2013) 418–442. [28] K. Mattsson, M.H. Carpenter, Stable and accurate interpolation operators for high-order multi-block finite-difference methods, SIAM J. Sci. Comput. 32 (4) (2010) 2298–2320. [29] K. Mattsson, F. Ham, G. Iaccarino, Stable and accurate wave propagation in discontinuous media, J. Comput. Phys. 227 (2008) 8753–8767. [30] K. Mattsson, J. Nordström, Summation by parts operators for finite difference approximations of second derivatives, J. Comput. Phys. 199 (2) (2004) 503–540. [31] K. Mattsson, F. Parisi, Stable and accurate second-order formulation of the shifted wave equation, Commun. Comput. Phys. 7 (2010) 103–137. [32] K. Mattsson, M. Svärd, M.H. Carpenter, J. Nordström, High-order accurate computations for unsteady aerodynamics, Comput. Fluids 36 (2006) 636–649. [33] K. Mattsson, M. Svärd, M. Shoeybi, Stable and accurate schemes for the compressible Navier–Stokes equations, J. Comput. Phys. 227 (4) (2008) 2293–2316. [34] Ken Mattsson, Summation by parts operators for finite difference approximations of second-derivatives with variable coefficients, J. Sci. Comput. 51 (2012) 650–682. [35] J. Nordström, Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation, J. Sci. Comput. 29 (2006) 375–404, http://dx.doi.org/10.1007/s10915-005-9013-4. [36] J. Nordström, M.H. Carpenter, Boundary and interface conditions for high-order finite-difference methods applied to the Euler and Navier–Stokes equations, J. Comput. Phys. 148 (1999) 341–365. [37] J. Nordström, M.H. Carpenter, High-order finite difference methods, multidimensional linear problems, and curvilinear coordinates, J. Comput. Phys. 173 (2001) 149–174. [38] J. Nordström, J. Gong, A stable hybrid method for hyperbolic problems, J. Comput. Phys. 212 (2006) 436–453. [39] J. Nordström, K. Mattsson, R.C. Swanson, Boundary conditions for a divergence free velocity–pressure formulation of the incompressible Navier–Stokes equations, J. Comput. Phys. 225 (2007) 874–890. [40] T.H. Pulliam, D.S. Chaussee, A diagonal form of an implicit approximate-factorization algorithm, J. Comput. Phys. 39 (1981) 341–363. [41] S. De Rango, D.W. Zingg, A high-order spatial discretization for turbulent aerodynamic computations, AIAA J. 39 (7) (2001) 1296–1304. [42] B. Sjögreen, High order centered difference methods for the compressible Navier–Stokes equations, Technical report 01.01, RIACS, NASA Ames Research Center, 2001. [43] B. Strand, Summation by parts for finite difference approximations for d/dx, J. Comput. Phys. 110 (1994) 47–67. [44] John C. Strikwerda, High-order-accurate schemes for incompressible viscous flow, Int. J. Numer. Methods Fluids 24 (1997) 715–734. [45] M. Svärd, On coordinate transformation for summation-by-parts operators, J. Sci. Comput. 20 (1) (2004). [46] M. Svärd, M.H. Carpenter, J. Nordström, A stable high-order finite difference scheme for the compressible Navier–Stokes equations, no-slip wall boundary conditions, J. Comput. Phys. 227 (May 2008) 4805–4824. [47] M. Svärd, K. Mattsson, J. Nordström, Steady-state computations using summation-by-parts operators, J. Sci. Comput. 24 (1) (July 2005) 79–95. [48] M. Svärd, J. Nordström, Stability of finite volume approximations for the Laplacian operator on quadrilateral and triangular grids, Appl. Numer. Math. 51 (1) (October 2004). [49] M. Svärd, J. Nordström, On the order of accuracy for difference approximations of initial–boundary value problems, J. Comput. Phys. 218 (October 2006) 333–352. [50] B. Yang, D. Gottlieb, J.S. Hesthaven, Spectral simulations of electromagnetic wave scattering, J. Comput. Phys. 134 (2) (1997) 216–230.