Energy stable numerical methods for hyperbolic partial differential equations using overlapping domain decomposition

Journal of Computational Physics 231 (2012) 5243–5265

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Energy stable numerical methods for hyperbolic partial differential equations using overlapping domain decomposition Adam Reichert a, Michael T. Heath a, Daniel J. Bodony b,⇑ a b

Department of Computer Science, University of Illinois at Urbana-Champaign, United States Department of Aerospace Engineering, University of Illinois at Urbana-Champaign, United States

a r t i c l e

i n f o

Article history: Received 15 April 2011 Received in revised form 16 February 2012 Accepted 13 March 2012 Available online 3 April 2012 Keywords: High order finite difference methods Overlapping domain decomposition Numerical stability Generalized summation-by-parts

a b s t r a c t Overlapping domain decomposition methods, otherwise known as overset grid or chimera methods, are useful for simplifying the discretization of partial differential equations in or around complex geometries. Though in wide use, such methods are prone to numerical instability unless numerical diffusion or some other form of regularization is used, especially for higher-order methods. To address this shortcoming, high-order, provably energy stable, overlapping domain decomposition methods are derived for hyperbolic initial boundary value problems. The overlap is treated by splitting the domain into pieces and using generalized summation-by-parts derivative operators and polynomial interpolation. New implicit and explicit operators are derived that do not require regularization for stability in the linear limit. Applications to linear and nonlinear problems in one and two dimensions are presented, where it is found the explicit operators are preferred to the implicit ones. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Domain decomposition methods offer an efficient way to solve initial-boundary value problems (IBVPs) for partial differential equations (PDEs) on general domains. In this approach, an IBVP is solved by splitting its domain into pieces, called subdomains, and then defining a new IBVP subproblem for each subdomain. The boundary conditions for each subproblem describe how that subproblem relates to its neighbors. An overlapping decomposition is one in which the subdomains have interior points in common, otherwise it is called non-overlapping. The idea has its origins with the work of Schwarz [1] who solved Laplace’s equation on a domain comprising an overlapping circle and square, like the keyhole domain in Fig. 1, as a sequence of problems that could be solved analytically. Domain decomposition has subsequently been used to transform an IBVP into a set of interrelated IBVPs that could be stored and solved on a computer. Early non-overlapping domain decomposition methods, called substructuring methods, were used to solve elliptic problems from structural analysis [2] by decomposing a large domain, whose size exceeded machine resources, into a sequence of smaller, interrelated ones. In one version of substructuring, the solution of one of these smaller linear systems gives the behavior of a physical substructure assuming its boundary is fixed, and the solution of another system gives the difference between the former solution and the state in which the boundaries of this substructure and its neighbors are allowed to move [3]. Comprehensive theoretical results are available for overlapping domain decomposition methods for elliptic problems [4–6]. For example, there is a provably convergent method for numerically solving

⇑ Corresponding author. E-mail addresses: [email protected] (A. Reichert), [email protected] (M.T. Heath), [email protected] (D.J. Bodony). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.03.003

5244

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

(a)

(b)

(c)

Fig. 1. (a) Keyhole domain X, with dark line around outside representing boundary @ X. (b) Overlapping domain decomposition comprising a square and a circle, with darker region representing overlap. (c) Non-overlapping decomposition.

Laplace’s equation on polygons that have been decomposed into overlapping pieces [7,8]. Others deal with elliptic equations on regions with curved boundaries [9] or consider finite element based approaches [10]. Overlapping domain decomposition methods have also been applied to time-dependent, nonlinear, fluid-dynamics problems (e.g., [11–16]), including aerodynamics and aeroacoustics problems (e.g., [17–23]) as a means to leverage high-order methods (usually finite difference) for use in complex geometries, and several tools are available to generate the inter-grid communication links and interpolants [24–38]. Adaptive overset methods have also been developed [39–43]. Even though overlapping methods are widely used, few proofs are available for their stability. Several conservative approaches have been developed [44–46] but without an explicit proof of stability and convergence. In some cases instabilities are observed when the primitive variables are interpolated between grids whereas the solution has improved robustness when the conservative variables are interpolated instead [46]. Some non-overlapping methods have been proven stable, including methods for the second-order wave equation [47,48], the advection equation in two spatial dimensions when one subdomain surrounds the other [49], and hyperbolic systems that place restrictions on the discretization of the subdomains [50,51]. Non-overlapping stability for the compressible Navier–Stokes equations on multiblock grids was determined by Nordström et al. [52]. Unfortunately, the non-overlapping proof techniques do not easily generalize to the overlapping case. Overlapping numerical methods have typically been proven stable on a method-by-method basis. For example, stability of one overlapping method is enforced by the addition of an empirically determined amount of artificial dissipation [53]. An overlapping domain decomposition method for the specific case of the advection equation in one spatial dimension using the Lax–Wendroff method has also been proven stable [54]. Another method for the advection equation has also been shown to converge [23]. While other authors have enforced stability of overlapping domain decomposition methods by the addition of artificial dissipation or proven that specific, low-order schemes are stable, we develop new, provably energy stable methods for time-dependent hyperbolic IBVPs on overlapping domains that are both high-order and include no artificial dissipation. We do this by using a class of easily constructible matrices called generalized summation-by-parts (GSBP) operators and showing that by linking portions of the operators with interpolating polynomials, strict energy stability can be proved. Several GSBP operators are constructed using a new, automated methodology that applies to both implicit and explicit operators. In Section 2 we briefly review necessary definitions before presenting existing energy stability methods and GSBP operators in Section 3. The new methodology is first developed in Section 4 in the context of the advection equation on overlapping domains before extending the method to systems of hyperbolic equations in Section 5. In Section 6 we present several numerical results for which proofs are not available. 2. Definitions and details 2.1. Hyperbolic systems of equations Consider the following notation for a k  k, real symmetric matrix A:  N A , the indicator of A, is a k  k matrix defined by

( fN A gi;j ¼

1 fAgi;j > 0; 0

fAgi;j 6 0:

 KA = diag(k1, k2, . . . , kk), where k1 P k2 P    P kk are the (necessarily real) eigenvalues of A.  Aþ ¼ TN KA KA T T and A = A  A+, where T is an orthogonal matrix that diagonalizes A,

A ¼ T KA T T :

ð1Þ

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

5245

The matrices A+ and A are uniquely defined: if T1 and T2 are two matrices that both satisfy (1), then the orthogonal matrix ðT T1 T 2 Þ commutes with ðN KA KA Þ. Let X  Rd and u : X  R ! Rk . Let A1, . . . , Ad be k  k real symmetric matrices. For a point z on oX, let nðzÞ 2 Rd be the outward pointing unit normal. Define B(z), the boundary matrix, as

BðzÞ ¼ n1 ðzÞA1 þ n2 ðzÞA2 þ    þ nd ðzÞAd : The system of partial differential equations under consideration is

ut ¼

d X Ai uxi

8 t > 0;

i¼1 þ

ðBðzÞÞ u ¼ f ðz; tÞ 8 t > 0; uðz; 0Þ ¼ gðzÞ

z 2 X; z 2 @ X;

ð2Þ

8 z 2 X;

where uxi ¼ @ xi u; f ðzÞ characterizes the boundary conditions, and g(z) gives the initial conditions. For k = 1, (2) describes the advection equation in d dimensions. We always assume that X, Ai, f, and g are such that the system is well-posed, as discussed in [55]. The IBVP (2) is called energy stable if the change in the norm of the solution can be bounded by an integral over the boundary,

kuk2 ¼

Z

uðz; tÞT uðz; tÞz; Z d kuk2 6 f ðz; tÞT CðzÞf ðz; tÞz; dt @X X

where C(z) is the pseudo-inverse of B+(z) defined by the Moore–Penrose conditions [56]. It follows that the existence of a solution implies the uniqueness of that solution. Moreover, the norm of the difference between any two solutions with different initial values never increases. We discretize the spatial operators @ xi in (2) using finite differences. The time evolution of v is then given by a system of linear ordinary differential equations obtained by the method of lines,

v_ ¼ Mv þ b; v ð0Þ ¼ g^;

ð3Þ

where v_ ¼ dv =dt and g^ are the initial data evaluated at the grid points. The matrix M represents a finite difference operator P that approximates the differential operator Ai @ xi and the boundary operator B+ in a form that will be made precise later. The inhomogeneous term b captures the effects of the boundary function f in (2). 2.2. Measures of error We use two measures of error. The truncation error is given by

^_  ðMu ^ þ bÞ; t¼u

ð4Þ

^ is the exact solution evaluated at the discrete grid points. Similarly, the cumulative error is defined by where u

^: e¼vu Our goal is to make e small in a sense described in the next section. 2.3. Convergence We will solve problem (2) by approximating the solution with a sequence of increasingly-refined overlapping composite b 1; X b 2 ; . . . be a sequence of progressivelygrids and will show that the cumulative error goes to zero in the natural norm. Let X i b refined grids discretizing X. Let ni be the number of points in X . Suppose that a numerical method produces a sequence of b ¼ fX b i gn forms a sequence of uniform grids, approximate solutions v1, v2, . . . , one for each discretization of the domain. If X i¼1 then we say the method converges to the true solution u at fixed time tf if

limkeðtf ÞkR ¼ 0; b X

ð5Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b as n ? 1. The norm kxk ¼ DVxT x is the natural where the limit is understood to be taken over the sequence of grids X R norm where DV is a measure of the volume of the grid cells. For example, for a one-dimensional uniformly spaced grid with spacing h, DV = h, while on a two-dimensional uniform grid DV = hxhy, the product of the one-dimensional grid spacings.

5246

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

2.4. Consistency The sequence of ODEs defined by M and b is consistent if

lim sup ktðtÞkR ¼ 0; b X t2½0;tf  b i , then the method is where t is the truncation error from Eq. (4). Further, if h maximum grid spacing for a discretization X consistent with order p if p

sup ktðtÞkR ¼ Oðh Þ

t2½0;tf 

as the grids are refined. 3. Energy stability, summation-by-parts, and simultaneous approximation term We say a scheme, or method, is energy stable in the H-norm if, for every system matrix M corresponding to the implementation of the scheme on a progressively-finer discretization, there is a norm H such that (a) H is a well-behaved norm in the sense of [57] and (b) M is negative semidefinite in the H-norm,

xT HMx 6 0 8x:

ð6Þ

Note that (6) is equivalent to all of the eigenvalues of the real symmetric matrix (HM + MTH) being non-positive and implies that the eigenvalues of M have non-positive real parts. An important theorem relates energy stability and convergence: Theorem 3.1. A consistent and energy stable method converges with bound

^ ðt f Þ  v ðtf ÞkR 6 ku

tf c1

! sup ktðtÞkH ;

ð7Þ

t2½0;t f 

where c1 is the constant from bound c1 kxk2R 6 kxk2H 6 c2 kxk2R ; 8x by norm equivalence. A proof can be found in [58]. Note that if M is consistent with order p, then

  u ^ ðtf Þ  v ðt f ÞR ¼ Oðhp Þ: We note that bound (7) is often not tight but making stronger assumptions on the structure of the product of H and M allows one to derive a slightly tighter bound [49]. For some specific cases, asymptotically tighter bounds can be derived. For example, [59] proves a tighter bound for the advection–diffusion equation with a particular approximation to the first and second derivatives. 3.1. Summation-by-parts As discussed in [58,60,61,57], a summation-by-parts (SBP) finite difference approximation to ox using the matrices {P, Q} satisfies the properties that (a) P is a well-behaved norm, (b) P1Q is an order p approximation to @ x, and (c) Q + QT = diag(1, 0, . . . , 0,1). When these conditions are satisfied the discrete version of integration-by-parts,

hP1 Qx; yiP ¼ xn yn  x1 y1  hx; P1 Q yiP

ð8Þ

holds. There are many ways to construct pairs of SBP matrices {P, Q}, and pairs have been constructed for which P1Q is up to a locally eighth-order approximation to @ x [58,60,61,57]. If P is diagonal, then the operator P1Q is termed explicit. Otherwise, it is termed implicit. It can be shown that assuming a certain structure for Q precludes P from being proportional to the identity [62]. A pair of matrices {Pc, Qc} satisfies the generalized summation-by-parts (GSBP) property [63–65] if it satisfies the first two e u ; 0; . . . ; 0; Q e l Þ where Q e u , and Q e l are symmetSBP properties (a) and (b) and a modified property, (c0 ) Q c þ Q Tc ¼ block diagð Q ric blocks. Clearly any pair of matrices satisfying the SBP property also satisfy the GSBP property. 3.2. Simultaneous approximation term The simultaneous approximation term (SAT) methodology is a process for enforcing boundary conditions such that the resulting method is energy stable. At an SAT boundary one or more penalty parameters s are introduced, which may not be directly related to the continuous problem, whose values are sought such that the resulting system is energy stable. The utility of the SAT methodology is that it reduces the problem of finding energy stable methods to enforcing that a small number of penalty parameters are set to the correct values. SAT methods have been derived for several families of equations

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

5247

[47,59,66,67]. For example, in [59] the authors construct an SAT method for the advection–diffusion equation using two penalty parameters. In the next section we present an SAT method constructed from GSBP matrices for a single subdomain and show how to extend it to an energy stable method for two overlapping domains. 4. Energy stable method for advection in one dimension In this and the next section we will solve scalar and vector advection problems on a one-dimensional domain X = (a, b) with a uniform discretization x = [x1, x1, . . . , xn]T with h = xi+1  xi for i = 1, . . . , n  1 the stepsize and x1 = a and xn = b. We start with the scalar case (k = 1 in (2)) and build to systems (k > 1). While one-dimensional methods are interesting in their own right, we analyze them in order to motivate their generalizations for nonlinear problems in higher dimensions. We follow the path set by others and prove that our methods are stable and convergent for linear problems and demonstrate empirically that they are stable for nonlinear ones. Consider the model advection problem on the interval domain X = (a, b),

ut þ ux ¼ 0;

uðx; 0Þ ¼ gðxÞ;

uða; tÞ ¼ f ðtÞ:

ð9Þ

This problem is well-posed [68]. 4.1. Energy stable single-domain method We find it convenient to term the triplet {P, Q, w} a GSBP system of approximation order p and size n such that (a) {P, Q} is a GSBP pair of dimension n  n and approximation order p, (b) w 2 Rn where the vector w is called the penalty vector, and (c) the matrix G = Q + wsT, where s ¼ ½1; 0; . . . ; 0T 2 Rn , is negative semidefinite. G is called the energy matrix. The entries of w will depend on penalty parameters and the GSBP system being used. A GSBP/SAT system enforces boundary conditions weakly and the penalty vector describes how an approximate solution is penalized for disparity between the exact and computed boundary conditions. For example, if {P, Q} is a pair of SBP matrices, w = ss, and s P 12, then {P, Q, w} form a GSBP system because, for any x,

xT Q x  sxT ssT x ¼

ð1  2sÞx21  x2n 6 0: 2

Using GSBP systems, an energy stable method for solving (9) on a uniform discretization can be stated succinctly: Theorem 4.1. Let {P, Q, w} be a GSBP system of approximation order p. Then method

v_ ¼ P1 Q v þ P1 wðv1 ðtÞ  f ðtÞÞ; v ð0Þ ¼ g^

ð10Þ

is energy stable in the P-norm and converges to the true solution of (9) at time tf with order p, p

keðt f ÞkR ¼ Oðh Þ: Let G be the energy matrix associated with {P, Q, w} and M be the system matrix for (10). The proof follows from the fact that PM = G. 4.2. Energy stable overlapping method Here we derive the corresponding energy stable method for the advection equation on an overlapped, one-dimensional domain. Again consider the model advection Eq. (9), but this time the domain X is decomposed into two parts: b 1 and X b 2 are uniform discretizations defined by X1 = (a = aL, bL) and X2 = (aR, bR = b). The sets X

xLi ¼ aL þ ði  1ÞhL

i ¼ 1; . . . ; n;

xRi

i ¼ 1; . . . ; m;

¼ aR þ ði  1ÞhR

where hL ¼ ðbL  aL Þ=ðn  1Þ; where hR ¼ ðbR  aR Þ=ðm  1Þ

ð11Þ

and illustrated in Fig. 2. The left grid, represented by the upper line, discretizes (a = aL, bL). The right grid, represented by the lower line, discretizes (aR, bR = b). Numbers above grid points give the ordering. When refining the grids, we move the right boundary of the left grid to keep the number of points in the overlap fixed; that is, bL = aR + (n + 1/2)hL, where n + 1 is the number of points in the overlap. Theorem 4.2 describes a new, energy stable, overlapping domain decomposition method that solves problem (9). For convenience, MATLAB-style notation is used to denote subsets of vectors,

fxg½i:j ¼ ½xi ; . . . ; xj T and we also use an unsubscripted norm kxk ¼

pffiffiffiffiffiffiffiffi xT x.

5248

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

1

2

3

4

n-3

n-2 n-1

aL

n bL

1

2

3

4

5

aR

m-1

m

bR

Fig. 2. Interval (a, b) discretized by two overlapping grids.

Theorem 4.2. Assume that X = (a, b) is discretized according to (11). Let {PL, QL, wL} be a GSBP system of approximation order p and size n, with energy matrix GL. Let {PR, QR, wR} be a GSBP system of approximation order q and size and m, with energy matrix GR. Let I 2 Rn be an interpolation vector satisfying r

^ 1 ðtÞ ¼ uðaR ; tÞ þ Oðh Þ: I Tu Consider the method 1 1 1 v_ 1 ¼ P1 v 1 ð0Þ ¼ g^1 ; L Q L v þ P L wL ðv 1  f ðtÞÞ; 1 1 T 1 2 2 2 v_ ¼ PR Q R v þ PR wR ðv1  I v Þ; v 2 ð0Þ ¼ g^2 :

ð12Þ

Assume the following e l Þ, where Q e l is sL  sL and sL P 1. Further, Q e l is positive definite.  Q L þ Q TL ¼ block diagð. . . ; 0; Q L L L e u ; 0; . . .Þ, where G e u is sR  sR and sR P 1. Further, G e u is negative definite.  GR þ GTR ¼ block diagð G R R R  The right penalty vector has a fixed number of entries of fixed size:

fwR gj ¼ 0;

j ¼ sR þ 1; . . . ; m;

 Rk ¼ C1; kw  R ¼ fwR g½1:sR  . where w  Most entries of I are zero, and the rest have bounded size:

I j ¼ 0; j ¼ 1; . . . ; n  sL ;   I  < C 2 ; where I ¼ fI g½nsL þ1:n .  Define 

a ¼

    e l kmax G eu kmin Q L R C 21 C 22

;

e l Þ and kmax ð G e u Þ denote the minimum eigenvalue of Q e l and the maximum eigenvalue of G e u , respectively. where kmin ð Q L R L R  Define s = min(p, q, r). Then method (12) is energy stable in the H-norm, where

H¼





PL

ð13Þ

aPR

and a is any scalar satisfying 0 < a 6 a⁄, and convergent in the sense of (5) with bounds s

ke1 ðtf ÞkRL ¼ Oðh Þ;

s

ke2 ðt f ÞkRR ¼ Oðh Þ:

ð14Þ

The parameter a is neither related to the original PDE (9), nor does it appear in method (12). This reflects that the norm matrix H is a theoretical tool used to show several properties. Its existence proves that method (12) is energy stable, and by Theorem 3.1, it follows that method (12) is also convergent. By the Lyapunov theorem, one can prove that all eigenvalues of the system matrix have non-positive real parts. We prove Theorem 4.2 by first proving an auxiliary lemma and then applying that lemma to method (12): Lemma 4.3. Let a, b, c, and b be real numbers. Assume a > 0 and c > 0. Define

fb ðx; yÞ ¼ ax2 þ bbxy  bcy2 : Then there exists b⁄ > 0 such that if 0 6 b 6 b⁄, then for all x; y 2 R; f b ðx; yÞ 6 0.

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

5249

The statement holds if b = 0 and b⁄ = 1. Assume b – 0. Observe that

fb ðx; yÞ ¼ ½x; yF

  x y

" ;

F¼

a

b 2b

b 2b

bc

# :

Let b⁄ = 4ac/b2. If 0 6 b 6 b⁄, then det F P 0. By Sylvester’s criterion, matrix F is negative semidefinite. Next, we apply Lemma 4.3 to method (12). We define the following vectors that represent the approximate solution on the rightmost nodes of the left domain and the leftmost nodes of the right domain,

l ¼ fv 1 g½nsL þ1:n ; r ¼ fv 2 g½1:sR  : Observe that

 R I T l: ðv 2 ÞT wR I T v 1 ¼ rT w If a 6 a⁄, then

    d T el 2 2 T T T eu el eu  kv k2H ¼ 2v T HMv 6 l Q L l þ 2ar wR I l þ ar G R r 6 kmin Q L klk þ 2akwR kkI kklkkrk þ akmax G R krk 6 0; dt where the last inequality follows from Lemma 4.3. This completes the proof of Theorem 4.2. A key aspect of Theorem 4.2 is that energy stability rests on being able to construct the a GSBP pair for the left (upwind) e l , of size sL  sL, is positive definite. The size of the interpolation vector I is tied to domain {PL, QL} whose lower right block Q L l e e l —so long as it is consistent and the size of the block Q L . The choice of the interpolation is arbitrary—and is not tied to Q L bounded. If the interpolation is based on polynomial interpolation, then any overlap that remains within the spatial support of I as the mesh is refined is acceptable. 4.3. Example We illustrate Theorem 4.2 by using a pair of GSBP matrices that are a locally fourth-order Padé approximation to the derivative on the interior of the domain [63–65] in conjunction with second-order interpolation. For convenience, they are reproduced in Appendix A. These operators were chosen because they have a simple structure with 2  2 upper and lower symmetric blocks. Assume that X has been discretized by two overlapping uniform grids, as in Fig. 2. Assume that the leftmost node of the right domain lies between the two rightmost nodes of the left domain,

xLn1 6 aR 6 xLn :

ð15Þ

Let h = 1  (aR  bL)/hL. For this choice of GSBP pair and second-order linear interpolation, we have

sL ¼ 2;

sR ¼ 2;

  5 1  R ¼ s ;  ; w 8 4 I ¼ ½h; 1  h;

" el ¼ Q L

1 4 1 4

1 4 3 4

# ;

 e u ¼ diag ð1  sÞ 5 ;  1 ; G R 4 4

where s is a penalty parameter. Therefore,

 R k2 ¼ s2 kw

 25 1 þ ; 64 16

kIk2 6 1;

   e u ¼ max ð1  sÞ 5 ;  1 ; kmax G R 4 4

  e l 0:146: kmin Q L

Letting s = 2, then a⁄ 0.0225. 4.4. Constructing GSBP pairs As mentioned in Theorem 4.2 the key to stability is being able to construct the GSBP pair of matrices {P, Q} where Q has e l is positive definite. Chertock [63] found two such operators, one of which is given in the property that its lower block Q L Appendix A, using an elaborate method that is difficult to generalize. (also Refs. [64,65].) In Appendix B we present an automated method for constructing a wider class of GSBP pairs {P, Q} and provide three examples. The first is an implicit, fourthorder approximation to @ x on the interior and third-order at the boundary. In this case, the upper symmetric block is 4  4, which admits the use of fourth-order interpolation. The second pair of matrices is implicit, sixth-order on the interior and

5250

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

fifth-order at the boundary. The lower symmetric block is 6  6 and admits sixth-order interpolation. We then construct an explicit, fourth order approximation to @ x on the interior that is second order at the boundary and generalizes the 2–4 SBP operator developed by Strand [61]. These pairs of matrices are constructed by simultaneously enforcing linear constraints ensuring that P1Q is a suitable first derivative operator and nonlinear constraints ensuring that P is positive definite and that the upper and lower symmetric blocks are indefinite and positive definite, respectively. We present the construction of the matrices, the final results, and a Fourier analysis of them in Appendix B. 5. Energy stable method for hyperbolic systems in one dimension In this section, we show how to extend the method presented in the previous section to an energy stable method for the hyperbolic system

ut ¼ Aux ; uðx; 0Þ ¼ gðxÞ;

ð16Þ

ðI  N A Þuða; tÞ ¼ fa ðtÞ; N A uðb; tÞ ¼ fb ðtÞ;

where A is a k  k matrix. We proceed by first introducing the reversing operator, a convenient tool for dealing with GSBP matrices. Next, we show how to extend the GSBP method from the previous section to scalar advection problems with arbitrary wave speeds. Finally, we solve the general hyperbolic problem (16). 5.1. Reverse of a matrix It is convenient to introduce the reversing operation for matrices. Let A be an m  n matrix. The matrix A], called the reverse of A, is the m  n matrix containing the elements of A permuted according to

fA] gi;j ¼ fAgmiþ1;njþ1 : For example,



1 2 3

]

 ¼

4 5 6

6 5 4 3 2 1

 :

The reversing operation satisfies the following properties, where A and B are any matrices:    

(A])] = A. (aA + bB)] = aA] + bB]. (AB)] = A]B]. If A is invertible, then A] is also invertible, with (A])1 = (A1)]. In this case we write A].

If A = A], then A is called reversible. If A = A], then A is called antireversible. Note that for any reversible square matrix A and any vector of the same dimension x,

xT Ax ¼ ðx] ÞT Aðx] Þ:

ð17Þ

The reversing operation behaves nicely with respect to discrete derivative operators, as the next theorem shows: Theorem 5.1. If D is an approximation to the derivative @ x of approximation order p, then D] is an approximation to the derivative @ x of the same approximation order. Let f(x) be a differentiable function. Define g(x) = f(x). Observe that

^f ¼ g^] ;

^f x ¼ g^] : x

Computing the truncation error associated with D] reveals that

D] ^f  ^f x ¼ D] g^] þ g^]x ¼ ðDg^  g^x Þ : ]

Because the natural norm R is reversible, it follows from (17) that

   ]^ ^  p D f  f x  ¼ kDg^  g^x kR ¼ Oðh Þ: R

We use the following corollary: Corollary 5.2. If {P, Q} is a pair of GSBP matrices of any approximation order, then P]Q] is an approximation to @ x of the same approximation order.

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

5251

5.2. General wave speed The first step in extending the energy stable method (12) to problem (16) is to consider the advection problem

ut  ux ¼ 0;

uðx; 0Þ ¼ gðxÞ;

uðb; tÞ ¼ f ðtÞ

ð18Þ

on the interval X = (a, b). The domain X is discretized according to (11). e describe an order r Theorem 5.3. Let {PL, QL, wL} and {PR, QR, wR} be GSBP systems of order p and q, respectively. Let I interpolation,

^ 2 ðtÞ ¼ uðbL ; tÞ þ OðhrR Þ: Ie T u Consider





] 1 ] ] T 2 1 v_ 1 ¼ P] v1 ð0Þ ¼ g^1 ; L Q L v þ P L wL v n  I v ;

 ] 2 ] ] 2 v_ 2 ¼ P] v2 ð0Þ ¼ g^2 : R Q R v þ P R wR v m  f ðtÞ ;

ð19Þ

e e ] , method (19) is energy stable in the H-norm, Under the same assumptions as Theorem 4.2 with I ¼ I where

"

e ¼ H

#

aP]L

ð20Þ

P]R

and a is any scalar satisfying 0 < a 6 a⁄, and converges to the solution of problem (18) with order min (p, q, r). Note that P]Q] approximates the derivative +@ x on both the left and right domains, which follows from Corollary 5.2. The proof follows from considering a left-moving wave as a right-moving wave under a change of variables. We now proceed to systems with arbitrary wave speed. Let {P, Q, w} be a GSBP system. For simplifying the notation we define the following:

D ¼ P1 Q  P1 wsT ; c ¼ P1 w; s ¼ ½1; 0; . . . ; 0T ; e ¼ P] Q ] þ P ] w] rT ; ~c ¼ P] w] ; r ¼ ½0; . . . ; 0; 1T : D

ð21Þ

The general advection equation in one dimension

ut  kux ¼ 0;

uðx; 0Þ ¼ gðxÞ;

ð1  N k Þuða; tÞ ¼ fa ðtÞ;

N k uðb; tÞ ¼ fb ðtÞ;

ð22Þ

where N k is the indicator of k defined earlier and ‘switches on’ the correct boundary condition. Let {PL, QL, wL} and {PR, QR, wR} e L; ~ be two GSBP systems. Using PL, QL, and wL, let DL ; cL ; D cL be the derivative operators and penalty vectors defined by (21) e R; ~ when P = PL, Q = QL, w = wL. The right operators are defined similarly. Using PR, QR, and wR, let DR ; cR ; D cR be the derivative e operators and penalty vectors defined by (21) when P = PR, Q = QR, w = wR. Further, assume fPL ; Q L ; wL g; fP R ; Q R ; wR g; I , and I meet all the criteria of Theorems 4.2 and 5.3. Let

"

  v_ ¼ kþ De þ k D v  fb ðtÞd~  fa ðtÞd; "

v ð0Þ ¼ g^;

~ ¼ D

eL D

~cL Ie

0m;n

eR D

#

D¼ "

;

~¼ d

0n ~cR

DL

0n;m

cR I

DR ;

"

# d¼

cL 0m

# ; ð23Þ

# :

Theorem 5.4. Let

e þ ð1  N þ ÞH Hk ¼ N kþ H k e defined by (13) and (20), respectively. Method (23) is energy stable in the Hk -norm and converges to the solution of with H and H (22) for any tf > 0. The key step in this proof is to perform a change of variables by rescaling time t0 = jkjt. The factor N kþ ‘‘switches on’’ the appropriate norm. 5.3. Systems of advection equations in one dimension A system of hyperbolic equations can be solved by decomposing it into characteristic variables and using the scalar method (23). We begin with a system that is already in characteristic form. Let A = K = diag(k1, k2, . . . , kk). Then

ut ¼ Kux ;

uðx; 0Þ ¼ gðxÞ;

ðI  N K Þuða; tÞ ¼ f ða; tÞ;

N K uðb; tÞ ¼ f ðb; tÞ

is a set of k decoupled advection problems. Let fa(t) = f(a, t) and fb(t) = f(b, t).

ð24Þ

5252

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

e and D and all penalty parameters be defined as in Theorem 5.4. Further, let Theorem 5.5. Let D

  ~ N þ þ ðH ðI  N þ ÞÞ H¼ H K K

ð25Þ

e defined by (13) and (20), respectively. Then method with H and H



v_ ¼ De Kþ v ð0Þ ¼ g^



þ ðD K Þ



v  ðd fa ðtÞÞ  ðd~ fb ðtÞÞ;

ð26Þ

is energy stable in the H-norm and converges to the solution of (24). Method (26) is energy stable and thus convergent because it is the sum of k energy stable, convergent, and independent components. Observe that v is a linear combination of k components, one for each wave speed ki. The time evolution of the component corresponding to ki is described by method (23) and is energy stable in the Hki -norm. The H-norm is the sum of the Hk1 ,. . ., Hkk norms, each applied to the corresponding component. Because the Hki -norm of each component is nonincreasing, their sum is also non-increasing. Because H is a symmetric permutation of a block diagonal matrix and each block e H is positive definite and is a well-behaved norm. Note also that the spectrum of the system matrix of (26) is is either H or H; the union of the spectra of the scalar systems (12) and (19). Because both are energy stable, all eigenvalues of the system matrix (26) have non-positive real parts. Next consider a general system (16) that is not in characteristic form. Let T be any orthogonal matrix that diagonalizes A,

T T AT ¼ KA : Let w(x, t) = TTu(x, t) be characteristic variables. Note that w solves

wt ¼ KA wx ; wðx; 0Þ ¼ T T gðxÞ; ðI  N KA Þwða; tÞ ¼ T T fa ðtÞ; N KA wðb; tÞ ¼ T T fb ðtÞ a system in the form of (24). Applying the diagonal method (26) to this new system and then transforming back gives the following result: e and D and all penalty parameters be defined as in Theorem 5.4. Further, let T be any orthogonal matrix that Theorem 5.6. Let D diagonalizes A and define







 e TN þ T T þ H T ðI  N þ ÞT T H¼ H K K e defined by (13) and (20), respectively. Then method with H and H



v_ ¼ D~ Aþ v ð0Þ ¼ g^



þ ð D A Þ



v  ðd fa ðtÞÞ  ðd~ fb ðtÞÞ;

ð27Þ

is energy stable in the H-norm and converges to the solution of (16). We have applied a similarity transformation to the system matrix from (26). Because similarity transformations preserve spectra, all eigenvalues of (27) also have non-positive real parts. Theorem 5.6 is the main result of this paper. It shows that for uniformly discretized overlapping domains an energy stable method utilizing GSBP pairs {P, Q} exists to a given order, provided the matrices {P, Q} can be constructed, as outlined in Appendix B . When the grids are non-uniform the results in [69,70] show that implicit SBP operators, where P is nondiagonal, fail to have a proof of energy stability. We expect these results to hold true for GSBP operators as well, so that on non-uniform meshes only the explicit GSBP operators would be (provably) energy stable. One comment is to be made regarding the extension of Theorem 5.6 to more than one dimension. As shown in [71], the generalization to multi-dimensional problems is straightforward when the problem is simultaneously diagonalizable, as in the simple two-dimensional scalar advection problem

@u @u @u þ kx þ ky ¼0 @t @x @y

ð28Þ

because the system’s characteristics all point in the same direction, (kx, ky), with increasing time. In this case the proof of Theorem 5.6 is identical so long as the coordinate ‘x’ is interpreted as being along the characteristic. When the system is not simultaneously diagonalizable, such as for the Euler equations in fluid mechanics or Maxwell’s equations in electromagnetics, our method of proof fails since there is a whole family of directions along which the characteristics propagate. We have not been able to prove energy stability for this class of problems, although our experience has been that they are stable in practice, as an example in the next section demonstrates for the two-dimensional Euler equations.

5253

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

6. Numerical results In this section, we present numerical results demonstrating the methods presented in this paper. For the one-dimensional problems the domain X is the open interval (1, 1) and the theoretical results are verified for scalar equations and for systems of equations. Examples for which we do not have theoretical results (non-linear, periodic domains, and twodimensions) are also given. In all cases we have chosen the overlap to consist of only one grid cell, regardless of the order of the method, for consistency amongst the examples. We stress that the overlap allowed is far more general and can lie anywhere within the spatial support of the interpolation vector I as defined first in Theorem 4.2. 6.1. Advection equation Consider the simple advection problem

ut ¼ ux

x 2 X;

uðx; 0Þ ¼ pðxÞ x 2 X;

ð29Þ

uð1; tÞ ¼ pð1  tÞ t > 0; where p(x) = sin(2px). The exact solution is u(x, t) = p(x  t). We test method (12) by solving the advection equation (29). The original domain X is decomposed into overlapping subdomains: X1 = (1, 0) and X2 = (hL/2, 1), where the distance of overlap diminishes with increasing resolution but the number of points in the overlap remains fixed. In order to achieve an energy stable method that converges with fourth-order on both domains, we must use at least fourth-order interpolation. By Theorem 4.2, energy stability requires the GSBP pair to have large lower symmetric blocks. We construct a new GSBP pair meeting this criterion in Appendix B.1. We have used this new GSBP pair in conjunction with fourth-order interpolation to achieve fourth-order convergence on both domains. The results are shown in Fig. 3(a). Appendix B.2 presents a new highorder GSBP pair with 6  6 lower symmetric blocks. We use these matrices and sixth-order polynomial interpolation to achieve sixth-order convergence on both domains. The results are shown in Fig. 3(b). 6.2. Periodic hyperbolic system The simple system

@u @u þ ¼ 0; @t @x @v @v  ¼0 @t @x

ð30Þ ð31Þ

with boundary conditions u(0, t) = av(0, t) and v(1, t) = bu(1, t) and initial conditions u(x, 0) = v(x, 0) = sin2px was used in [57] as a ‘severe’ test of stability because, when jaj = jbj = 1, the L2 energy of the system is constant for all time. We consider a twogrid version of this problem by splitting the problem into two domains x1 2 [1, 0] and x2 2 [hL/2, 1]. Observe that the size of the overlap region reduces as hL reduces. On the boundaries at x = 1 and x = 1 the usual SAT boundary conditions are applied while the overlap uses interpolation and the special block structure of Q for the left- and right-moving waves. Theoretical stability results are not available for this problem since the periodicity invalidates the method of proof used in Theorems 4.2 and 5.3. Instead, we report on the eigenvalues of the system matrix M, since these determine long time

10−4

.

10−8

right

10−5

10−7 right left

left

10−9

10−6

10−10

4

10−7 1 10−8

10−3

6 1

10−2

(a)

10−1

10−11 −3 10

10−2

10−1

(b)

Fig. 3. (a) Local convergence plot for problem (29) using method (12) with fourth-order Padé GSBP pairs, showing errors ke1 kPL and ke2 kPR as functions of mesh spacings hL and hR, respectively. (b) Same, using sixth-order Padé GSBP pairs.

5254

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

3 1

(a)

(b)

Fig. 4. (a) Convergence rate of explicit GSBP operator (h, left domain; ⁄, right domain) and (b) spectrum of eigenvalues of system matrix.

stability and, by the Lyapunov theorem, indicate that if the system is Lyapunov stable then there exists a norm H in which the system is energy stable. The implicit GSBP4 operator is unstable for this problem because of the coupling between the x = 1 and 1 boundaries, periodicity, and the overlap, while the explicit GSBP operator in Appendix B.3 is stable. Fig. 4 shows the convergence rate and eigenvalue spectrum for the stable case. (The right-most eigenvalue in Fig. 4(b) remains in the lefthalf-plane as m, n ? 1 so that the method remains stable as the grid is refined.) Thus, while the implicit GSBP operators are stable for overlap problems with inflow and outflow boundary conditions, they may not be stable when there is some periodicity in the system. The explicit operators remain stable, however, because the norms P do not couple all of the boundaries together. 6.3. Inviscid Burgers’ equation Nonlinear problems of the form (16), where A = A(u), can be solved using either explicit or implicit GSBP pairs. However, the use of explicit GSBP pairs creates an opportunity for efficiency. For many explicit GSBP pairs {P, Q}, the row of the discrete derivative operator P1Q that approximates the derivative at the interior point xi and the corresponding row of P]Q] are the same. Such an operator allows a discrete approximation to Aux jxi to be calculated without diagonalizing A. We have constructed such an explicit GSBP pair, which is shown in Appendix B.3. The operator P1Q is a fourth-order approximation to @ x in the interior and second-order at the boundary, and we have applied this globally third order explicit discrete derivative operator to the inviscid Burgers’ equation

1 0.8 right

0.6 left

0.4

3 1

0.2 3

0 −1

1

−0.5

0

(a)

0.5

1

(b)

Fig. 5. (a) Approximate solution to problem (32) at time t = 0.2 (b) local convergence plot for problem (32) using nonlinear extension of method (12), showing errors ke1 kPL and ke2 kPR as functions of mesh spacings hL and hR, respectively.

5255

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

0.6 1 0.8

0.5

0.6 0.4 0.2

0.4

0 −1 −0.8 −0.6 −0.4 −0.2

0

0.2 0.4 0.6 0.8

1

−0.1

0

(a)

0.1

0.2

(b)

Fig. 6. Contours of p  p1 showing (a) convection of vortex over time and (b) as it crosses the overlap interface.

ut ¼ uux ;

ð32Þ

uðx; 0Þ ¼ expð10x2 Þ;

where the exact solution to the Cauchy problem provides the boundary conditions. The results are shown in Fig. 5. As shown in Fig. 5(b), we achieve third order convergence on both the left and right subdomains, as expected.

6 coarse

4 2

1

fine

0

0

−2 med.

−4 −6 −8

−1 0

0.5

1

1.5

0

0.5

(a)

1

1.5

(b) −2

0.015

10

0.010

−3

10 0.005

−4

0

10

3.9

−0.005 −5

10

−0.010 −0.015

1

−6

0

0.5

1

(c)

1.5

10

−3

10

−3

10

−3

10

(d)

Fig. 7. (a) Scaled error (pnum  pexact)/dp vs. time at (1/2, 1/2). (b) Magnification of (a). (c) Scaled error vs. time at (1/2, 1/2). (d) Scaled error vs. mesh spacing at (1/2, 1/2) and peak reflection time t 0.8716.

5256

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

6.4. Euler equations We apply the method presented in this paper to a two-dimensional vortex governed by the non-linear Euler equations. In the reference frame of the vortex, fluid rotates around a stationary point where the pressure is locally a minimum. In the reference frame used for this computation, the vortex convects with a mean flow at speed v1 to the right. Distances are normalized by a characteristic length L so that the computational domain is X = [1, 1]  [0, 1], velocities are normalized by the speed of the vortex v1, and time is normalized by the characteristic time L/v1. The free-stream pressure p1 and free-stream density q1 are set to 1, and we use a ratio of specific heats c = 1.4. With these choices, the vortex is moving with unit speed, pffiffiffi has Mach-number 1= c 0:845, and starts at point x = (3/4, 1/2) at time t = 0. Initial and boundary data are provided by the analytical solution found in [72]. We discretize the left and right subdomains with overlapping, uniform composite grids with node counts given by the following table: Left

Coarse Medium Fine

Right

nx

ny

nx

ny

43 83 163

45 85 165

86 166 326

45 85 165

The leftmost edge of the right discretization always lies halfway between the rightmost two columns of the left discretization. The fourth-order GSBP pair presented in Appendix B.1 is used to approximate the derivatives @ x and @ y on both the left and right domains. Fourth-order polynomial interpolation is used at the interfaces of the two regions. Fig. 6 presents the results of the third simulation by plotting the pressure contours both before and after the vortex has crossed the interfaces, which are represented by the hashed lines. Observe that circular contour lines are retained. We next quantify the error in our method by measuring the errors at points away from boundaries and interfaces. These three discretizations were chosen because each has nodes placed at cL = (1/2, 1/2) and cR = (1/2, 1/2). We investigate the error in the computed pressure, ðpnum  pexact Þ=dp, where pnum is the computed pressure, pexact is the pressure given by the analytical solution, and 1  dp is pressure at the vortex center. Fig. 7(a)–(c) show the scaled errors at cL and cR as functions of time for each discretization. In Fig. 7(a), the increase in the magnitude of the error near time t = 1/4 corresponds to discretization error as the vortex moves over cL. Similarly in Fig. 7(c), the larger error near t = 5/4 corresponds to the discretization error as the vortex moves over cR. All simulations show a peak in the scaled error at cL and t 0.8716. Measuring this peak 3:9 shows the error reflected back into the domain from the interface. Fig. 7(d) shows that this error is Oðhx Þ, where hx is the xspacing between adjacent nodes in the left domain. 7. Conclusions Overlapping domain decomposition methods are applicable to a wide range of problems but are prone to instability unless numerical diffusion or some other form of regularization is used. To address this, we derived high-order, energy stable, overlapping domain decomposition methods for hyperbolic initial-boundary value problems that do not require artificial dissipation. First, we solved a model advection equation by splitting its domain into pieces and then applying new generalized summation-by-parts derivative operators and polynomial interpolation. We proved that this method is energy stable and convergent by first proving a lemma stating that a class of functions is negative semidefinite and then using this lemma to show that our method is energy stable in an infinite class of norms. Theorems 4.2 and 5.6, the main results of this paper, showed that high-order methods can be constructed by using high-order GSBP derivative operators and stated the conditions that these operators must satisfy. Appendices B.1 and B.2 present new GSBP derivative operators meeting these criteria that are fourth-order and sixth-order implicit methods, along with a fourth-order explicit method. A novel scheme for automatically constructing GSBP operators was also presented that extends the method originally developed in [63–65]. We applied these methods to several numerical examples to demonstrate the theoretical results. We also applied the new methods to problems for which a proof is not available. For example, both inviscid Burgers’ equation in one dimension and the Euler equations in two dimensions demonstrated empirical convergence. We observed that explicit GSBP derivative operators create an opportunity for efficiency by reducing the number of required eigenvalue decompositions and derived such an operator in Appendix B.3 that can be applied to non-uniformly stretched grids. Generalized summation-by-parts operators can be used to construct high-order, energy stable overlapping domain decomposition methods for both linear and nonlinear problems that do not require artificial dissipation. For problems of practical interest, which are likely non-linear and possibly periodic, the implicit (block) norm P GSBP operators are not energy stable (and may not be Lyaponuv stable); however, the diagonal norm P operators continue to be useful and are, therefore, preferred.

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

5257

Acknowledgments This work was supported, in part, by the US Department of Energy through the University of California under subcontract B523819, and by NASA (Subsonic fixed wing NRA, Dr. L. Hultgren, monitor). Appendix A. GSBP4 GSBP4 is a pair of GSBP matrices that are a locally fourth-order approximation to the first derivative on the interior and third-order at the boundary [63–65]. Note the factor of h multiplying P:

2

79

11 36 13 9 35 288 5 144

6 288 6 11 6 36 6  25 6 288 6 5 6 288 6 6 6 6 6 P ¼ h6 6 6 6 6 6 6 6 6 6 6 4 2

25  288 35 288 287 288 1 4

3

5 288 5 144 1 4 1 4

1 1 4

1 1 4 .. .. .. . . . 1 4

1

1 4

1 4

1 1 4 1 48 1 96

 58

451 576 1 8  181 192 5 144

6 595 6  576 6 6 29 6 144 6 25 6  576 6 6 6 6 6 Q ¼6 6 6 6 6 6 6 6 6 6 6 4

29  144

1 48 53 288 10 9 7 36

47 288

3

0

25 576 5  144 427 576

 427 576

0

3 4

 34

0 .. .

181 192

1 4 287 288 53 288 13  288

7 7 7 7 7 7 7 7 7 7 7 7 7; 7 7 7 7 1 7 96 7 13 7 7  288 7 7 7 36 5

3 4

..

.

..

 34

0

3 4

 34

0

143 192

1  48

 143 192 1 48 5  192

0

163 192 1 8  29 64

.

 163 192 5 48

7 7 7 7 7 7 7 7 7 7 7 7 7: 7 7 7 7 5 7 192 7 5 7 7  48 7 45 7 64 5 3 8

e u ; 0; . . . ; 0; Q e l Þ, where Explicit calculation reveals that Q þ Q T ¼ block diagð Q

" eu ¼ Q

 14 "

el ¼ Q

 54  14

1 4 1 4

1 4 3 4

1 4

# ;

# :

Define the penalty vector

 T 5 1 w ¼ s ;  ; 0; . . . ; 0 : 8 4 If s P 1, then {P, Q, w} is a GSBP system. All conditions in this paper are satisfied if s > 1. Appendix B. Matrices for overlapping domain decomposition B.1. Low-Order Implicit GSBP We present one possible procedure for constructing a pair of GSBP matrices {P, Q} that is a locally third-order approximation to @ x at the boundaries and fourth-order on the interior. The algorithm has three high-level steps:

5258

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

(a) Choose a Padé central difference approximation. (b) Select sizes of left and right boundaries. (c) Construct boundary blocks of P and Q. The problem of determining a GSBP pair for which P1Q approximates the derivatives with a given order of accuracy does not have a unique solution. We construct a solution by restricting ourselves to GSBP pairs in which both matrices are banded with corrections to the upper and lower diagonal blocks. The coefficients of the Padé scheme chosen in (a) determine the entries in the bands, and the sizes of the boundary blocks selected in (b) determine the sizes of the corrections. This three-step procedure is similar to that followed by other authors [58,61,57,63–65], but we complete step (c) in a new way. Whereas other authors analytically solve large systems, we use numerical optimization to improve an initial guess iteratively, in conjunction with a perturbation step that ensures that the pair {P, Q} satisfy in exact arithmetic a set of consistency constraints, which are defined later. The first step is to select an appropriate Padé approximation. We desire P1Q to be a fourth-order approximation on the interior. In order to satisfy this requirement, we select the fourth-order Padé approximation

h 0 h 3 3 0 f ðxiþ1 Þ þ hf ðxi Þ þ f 0 ðxi1 Þ f ðxiþ1 Þ  f ðxi1 Þ; 4 4 4 4

ðB:1Þ

where h = xi+1  xi. For simplicity, we have used the same Padé approximation used in the construction of GSBP4, but we emphasize that this is just one of many possible fourth-order choices [73]. Next, the sizes of the boundary blocks must be chosen. We choose all boundary blocks to be 4  4. Again, one could use differently sized blocks. The choices of Padé approximation and sizes of the boundary blocks assign the following structures to P and Q,

2

p1 6 6 p2 6 6 6 p3 6 6 6 p4 6 6 6 6 6 6 P ¼ h6 6 6 6 6 6 6 6 6 6 6 6 6 4 2

q1

6 6 q5 6 6 6 q9 6 6 6 q13 6 6 6 6 6 6 Q ¼6 6 6 6 6 6 6 6 6 6 6 6 6 4

3

p2

p3

p4

p5

p6

p7

p6

p8

p9

p7

p9

p10

1 4

1 4

1

1 4

..

..

.

1

1 4

1 4

p11

p12

p13

p12

p15

p16

p13

p16

p18

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7; 7 7 7 7 7 p14 7 7 7 p17 7 7 7 p19 7 5

p14

p17

p19

p20

..

.

1 4

.

3

q2

q3

q4

q6

q7

q8

q10

q11

q12

q14

q15

q16

3 4

 34

0

3 4

..

..

.

 34 q17

q18

q19

q21

q22

q23

q25

q26

q27

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7; 7 7 7 7 7 q20 7 7 7 q24 7 7 7 q28 7 5

q29

q30

q31

q32

.

..

 34

0

. 3 4

where the pi are unknowns in P and qi in Q. Observe that the coefficients of the Padé approximation (B.1) have become the values in the bands of P and Q. Also note that the sizes of the boundary blocks determine the number of unknowns. We now solve for these unknowns. Define

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

2

p11

6p 12 b¼6 P 6 4 p13 p14 1

P

p14

3

2

p12

p13

p15

p16

p16

p18

p17 7 7 7; p19 5

p17

p19

p20

q17

6q 21 b ¼6 Q 6 4 q25 q29

q20

5259

3

q18

q19

q22

q23

q26

q27

q24 7 7 7: q28 5

q30

q31

q32

Q being a third-order approximation to the derivative at the right boundary is equivalent to enforcing that

3 2 2 3 1 0 0 0 1 2 3 0 1 0 0 7 6 6 7 6 61 1 1 60 1 2 3 7 160 0 0 07 7 b6 6 7 b6 P 7 ¼ Q6 7þ 6 6 61 2 4 6 0 1 4 12 7 4 6 0 0 0 0 7 5 4 4 5 4 2

0 0

0 1 6 27

0

1 3 9 27

0

3

3 2 1 1 1 1 7 7 6 7 6 1 7 7 360 0 0 0 7 7þ 6 7; 460 0 0 0 7 8 7 5 5 4 0

0

0

0

0

which can be rewritten,

2

 11 6

3

 32

 12

1

3

1 3

2 23 12

7 6 6  16 7 7 60 7þ6 6 1 7 7 60 1 12 3 5 4 3 11 3 0 2 6

6 6 1 6 b ¼P b6 3 Q 6 1 6 6 4  13

 17 8

5 4

0

0

0

0

0

0

7  24

3

7 0 7 7 7: 7 0 7 5 0

ðB:2Þ

b. b in Eq. (B.2) uniquely determines Q Note that the choice of P b þQ b T Þ are symmetric positive definite. b such that both P b and ð Q So far, we have proceeded as in [63,64]. Next we obtain a P Work by other authors places more restrictions on system (B.2), uses symbolic math software to determine the number of b satisfying the desired propb and Q degrees of freedom, and then gives choices for these degrees of freedom that result in P b 0 . To do this, we create the formal erties. Instead, we use a numerical optimizer to improve iteratively a random initial guess P maximization problem

b trial ¼ arg max½min kðD þ DT Þ; P

ðB:3Þ

s:p:d:C

b ! D and P b ! C in Eq. (B.2), and k(X) is where C is a 4  4 real symmetric, positive definite matrix, D is defined by letting Q b 0 , we obtain the solution the spectrum of X. After solving (B.3) with MATLAB’s optimizer using the random initial guess P

2

0:7181

0:1581

0:7142

0:2963

3

7 6 6 0:1581 4:4364 5:2697 1:8557 7 7 b trial 6 P 7: 6 6 0:7142 5:2697 9:3293 3:7795 7 5 4 0:2963

1:8557

3:7795

2:623

b are defined to be rational approximations of the corresponding elements of P b trial , The elements of P

2

163 227

46  291

927 1298

8  27

3

7 6 6  46 244 180 7  469 6 291 55 89 97 7 b 7: 6 P¼6 927 793 257 7 7 6 1298  469  89 85 68 5 4 8 180 257 160  27  68 97 61 b is calculated according to Eq. (B.2), yielding Matrix Q

2

1007857804 1157516811

6 6  22954787 6 b ¼ 6 8546670 Q 6 8735503 6 2678010 4  102958735 65181672

 1080757783 1028903832

324312701 257225958

 3066198755 9260134488

2545369 474815

 3350641 949630

3646957 4273335

 401461839 39277480

47437449 4909685

 318475391 117832440

21354653 3621204

 54017093 7242408

102364409 32590836

b þQ b T Þ are positive definite, b and ð Q Note that P

b f0:500; 0:781; 1:578; 14:248g; kð PÞ b þQ b T Þ f0:258; 0:258; 1:784; 35:7689g: kð Q

3 7 7 7 7: 7 7 5

5260

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

It is not necessary perform an analogous procedure for the upper blocks. Because the upper blocks of GSBP4 satisfy the necessary requirements, we assign the upper blocks of P and Q to be the corresponding values of the GSBP4 matrices. The final matrices are

2

79 288 6 11 6 36 6  25 6 288 6 5 6 288 6

6 6 6 P ¼ h6 6 6 6 6 6 6 6 6 6 4

11 36 13 9 35 288 5 144

25  288 35 288 287 288 1 4

3

5 288 5 144 1 4

1 1 4

1 4 1 1 4 .. .. .. . . . 1 1 4 1 4

1 4 163 227 46  291 927 1298 8  27

46  291 244 55  469 89 180 97

927 1298  469 89 793 85  257 68

7 7 7 7 7 7 7 7 7 7 7; 7 7 7 7 8 7  27 7 180 7 7 97 7  257 5 68 160 61

3

2

451 29 25  58  144 576 576 1 181 5 6  595  6 576 8 192 144 6 29 427 0 6 144  181 192 576 6 25 5 427 3 6 0 4 6 576 144  576 6 3 3  0 6 4 4 6 6 . . .. .. ... Q ¼6 6 6  34 0 6 6 6  34 6 6 6 6 4

3 4 1007857804 1157516811  22954787 8546670 8735503 2678010  102958735 65181672

 1080757783 1028903832 2545369 474815  401461839 39277480 21354653 3621204

324312701 257225958  3350641 949630 47437449 4909685  54017093 7242408

7 7 7 7 7 7 7 7 7 7 7 7: 7 7 7 7 3066198755 7  9260134488 7 7 3646957 7 4273335 7 318475391 5  117832440 102364409 32590836

Define the penalty vector

 T 5 1 w ¼ s ;  ; 0; . . . ; 0 : 8 4 If s P 1, then {P, Q, w} is a GSBP system. All conditions in this paper are satisfied if s > 1. We next present a Fourier analysis of the discrete derivative operator P1Q defined by the above matrices. This type of analysis provides an effective way to quantify the resolving characteristics of a finite difference operator. Introductions to this topic as well as further reading can be found in [60,74]. We follow [73] in the details of our analysis. A uniform discretization of the domain X = [0, 2p) is the set n equally spaced points

2p ; n xi ¼ ði  1Þh;

h¼

i ¼ 1; . . . ; n:

ðB:4Þ

Let f : X ! R. The vector ^f , called the injection of f, is defined by

^f ¼ ½f ðx1 Þ; . . . ; f ðxn ÞT : The Fourier basis is the collection of functions g k ðxÞk 2 Z, where

g k ðxÞ ¼ expðikxÞ and i ¼

pffiffiffiffiffiffiffi 1. The elements of ^f can be written in terms of the Fourier basis according to

^f ¼ j

U X ~f g ðx Þ; k k j k¼L

where, if n is odd then L ¼  ðn1Þ and U =  L, and if n is even then L ¼  2n þ 1 and U = L + 1. The wavenumber w 2 [p, p] is 2 a scaled index defined by w ¼ 2pn k. For convenience, we define gw = gk(w). Note that

5261

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

 @g w  w ¼ i g w ðxj Þ: h @x xj An approximate first derivative D 2 Rnn defines a modified wavenumber w0j ðwÞ by

fDg^w gj ¼ i

w0j h

g w ðxj Þ:

The exact differentiation operator @ x corresponds to w0j ðwÞ ¼ w for all j. The value Reðw0 Þ  w measures dispersive error, and Imðw0 Þ measures dissipative error [73]. Fig. B.8 plots the modified wavenumber w0 vs. the wavenumber w for the discrete derivative operator P1Q at the boundary point x49 and the near-boundary point x39. Fig. B.8(a) shows that the GSBP modifications introduce dissipative error, while Fig. B.8(b) shows that this error has diminished in the interior points. In the next sections, we present two new pairs of GSBP matrices that can also be used to create energy stable methods for overlapping subdomains. The same procedure is followed except that different Padé approximations and sizes of boundary blocks are used. B.2. High-order implicit GSBP The following pair of GSBP matrices can be used to construct an approximation to the first derivative that is locally sixthorder on the interior, fifth-order at the boundaries:

2

P¼

6965509 27993600 6 30563 6 77760 6 6  281 6 1296 6 51007 6 349920 6 h6  19075 6 373248 6 19273 6 2332800 6

30563 77760 10875041 5598720 4159  15552 171491 466560 209207  1399680 45253 1866240

6 6 6 6 6

6 6. 6 .. ... ... 6 6 1 6 1 3 6 6 1 6 3 6 6 6 6 6 6 6 6 6 4 2

 23 6  7412017 6 5598720 6 1708109 6 2799360 6 6  38413 6 93312 Q ¼6 418787 6 2799360 6 6  27143 6 1119744 6 6 6 6 6 6 6. 6 .. ... 6 6 1 6  36 6 6 6 6 6 6 6 6 6 6 6 6 4

281  1296 4159  15552 1786169 1399680 10219 46656 27791 466560 7471  699840

1 3 298 321 34 91 2 47 1 377 13 259 2 67

5545777 5598720 1 6 246851  174960 999373 2799360 230747  1866240 11383 559872

34 91 28 43 29 90 3 37 9  103 1 1373

1708109  2799360

51007 349920 171491 466560 10219 46656 1382789 1399680 1 3

2 47 29 90 83 199 5 28 19 84 4 135

0

1009613  1399680 0 252707  2799360  129581 174960 5071 466560

27791 466560 1 3 5603021 5598720 1 3

1 377 3 37 5 28 87 163 1  46 1  209

.

..

 79

0

1  36

 79 1  36

.

13 259 9  103 19 84 1  46 3 5 4 51

418787  2799360 230747 1866240 252707 2799360 129581 174960

0

98947  2799360  4351153 5598720 1  36

..

..

19273 2332800 45253 1866240 7471  699840

209207  1399680

0

38413 93312 999373  2799360 1009613 1399680

246851 174960

19075  373248

0 1 3 27992569 27993600 1 3

7 7 7 7 7 7 7 7 2 67 7 7 1 7 1373 7; 7 4 135 7 7 1 7  209 7 7 4 51 5

3 7 7 7 7 7 7 7 7 7 7 7 1 7 3 7 7 1 1 7 3 7 .. .. .. 7 . . .7

4 45

27143 1119744 11383  559872 5071  466560 98947 2799360 4351153 5598720

0  79 .. .

3 7 7 7 7 7 7 7 7 1 7 36 7 7 7 1 7 9 36 7 7 1 7 0 9 36 7 .. .. .. .. 7 7 . . . .

.

7 1 9 36 160221934720031 1823322782801879 5828468439036149 1039568760741888 4851320883462144 14553962650386432 585283142126  29770625845139  1595135592781 36854611630200 921365290755 5528191744530 3997079 47879359 6403651  28283472  141417360 26515755 574092279179 53542596941 137761097569  18361750527120  1836175052712  306029175452 885990713 39157283 88129963  12516848400  625842420 312921210 17764323418 9593473514 49077364697  220619715525  132371829315 44123943105

1274863239037043 278292382087255  6823835967946752  232863402406182912 2599134554897  1633609455661 3685461163020 7370922326040 282899  35354340 278915914831 918087526356  209162939 625842420 18284546732 44123943105

19249987 141417360 904545390139 3672350105424 246728177 2503369680  17929655363 44123943105

7 7 7 7 7 7 7 1531531867889623 7 24256604417310720 7 7 145483789943 7 27640958722650 7: 7 46780673 7 424252080 7 120976857431 7  3060291754520 7 7 1007979137 5 3129212100 149454320494 661859146575

5262

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

3

3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0 −0.5 0

0.5

1

1.5

2

2.5

3

0

0

0.5

1

(a)

1.5

2

2.5

3

(b)

Fig. B.8. (a) Modified wave number vs. wave number at boundary point x49. (b) Same at near-boundary point x39.

10

3

5

2.5

0

2

−5

1.5

−10

1 0.5

−15 0

0.5

1

1.5

2

2.5

3

0 0

0.5

1

(a)

1.5

2

2.5

3

(b)

Fig. B.9. (a) Modified wave number vs. wave number at boundary point x49. (b) Same at near-boundary point x39.

3

3

2.5

2.5

2

2

1.5

1.5

1 0.5

1

0

0.5

−0.5 0

0.5

1

1.5

(a)

2

2.5

3

0 0

0.5

1

1.5

2

2.5

3

(b)

Fig. B.10. (a) Modified wave number vs. wave number at boundary point x49. (b) Same at near-boundary point x39.

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

5263

Define the penalty vector

 T 2 1 w ¼ s ;  ; 0; . . . ; 0 : 3 3 If s P 1, then {P, Q, w} is a GSBP system. All conditions in this paper are satisfied if s > 1. Fig. B.9 plots the modified wavenumber vs. the wavenumber for the discrete derivative operator defined by P1Q. Fig. B.9(a) shows that this operator introduces a dissipative error near the boundary, while Fig. B.9(b) shows that the dissipative error is much smaller in the interior. B.3. Explicit GSBP This final pair of GSBP matrices can be used to construct an approximation to the first derivative that is locally fourthorder on the interior, second-order at the boundaries:

P ¼ h diag 2

 12

 17 59 43 49 7 1 1 ; ; ; ; 1; . . . ; 1; ; ; ; 48 48 48 48 8 2 8 59 96

3

1 1  12  32

6 59 59 0 6  96 96 6 1 59 6 0 6 12  96 6 1 0  59 6 32 96 6 1 6 6 12 Q ¼6 6 6 6 6 6 6 6 6 6 4

0 59 96

0

1  12 2 3

 23 .. .

0 .. . 1 12

1  12 2 3

..

.

1  12 .. .

 23

0

1 12

 23 1 12

..

.

2 3 3 16  12 1 16

1  12 5 12 1 4  14

7 7 7 7 7 7 7 7 7 7 7: 7 7 7 7 7 7 1 7  48 7 1 7 6 5 3 16

Define the penalty vector

 T 1 w ¼ s ; 0; . . . ; 0 : 2 If s P 1, then {P, Q, w} is a GSBP system. All conditions in this paper are satisfied if s > 1. Fig. B.10 plots the modified wavenumber vs. the wavenumber for the discrete derivative operator defined by P1Q. Fig. B.10(a) shows that this operator introduces a dissipative error near the boundary. The derivative operator P1Q approximates the first derivative @ x in the interior using a fourth-order centered finite-difference approximation. Such approximations do not introduce dissipative error at the interior points. Fig. B.10(b) confirms this. Note that Imðw039 ðwÞÞ ¼ 0 for all w. References [1] H.A. Schwarz, Über einen grenzübergang durch alternierendes verfahren, Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich 15 (1870) 272–286. [2] B.F. Smith, P. Bjørstad, W. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, 2004. [3] J.S. Przemieniecki, Theory of Matrix Structural Analysis, Dover, 1985. [4] T.F. Chan, T.P. Mathew, Domain decomposition algorithms, Acta Numerica 3 (1994) 61–143. [5] A. Quarteroni, A. Valli, Domain Decomposition Methods for Partial Differential Equations, Oxford University Press, 1999. [6] A. Toselli, O.B. Widlund, Domain Decomposition Methods – Algorithms and Theory, Springer, 2005. [7] E.A. Volkov, A mesh method for finite and infinite polygons and error bounds in terms of known quantities, Proceedings of the Steklov Institute of Mathematics 96 (1970) 187–234. [8] E.A. Volkov, The method of composite meshes for finite and infinite regions with piecewise smooth boundaries, Proceedings of the Steklov Institute of Mathematics 96 (1970) 145–185. [9] G. Starius, Composite mesh difference methods for elliptic boundary value problems, Numerische Mathematik 28 (1977) 243–258. [10] P. Le Tallec, T. Sassi, M. Vidrascu, Three-dimensional domain decomposition methods with nonmatching grids and unstructured coarse solvers, in: C. Huneke (Ed.), Domain Decomposition Methods in Scientific and Engineering Computing, vol. 180, AMS, 1994, p. 61. [11] R. Magnus, H. Yoshihara, Inviscid transonic flow over airfoils, AIAA Journal 8 (1970) 2157–2162. [12] M.M. Rai, Navier–Stokes simulations of rotor/stator interaction using patched and overlaid grids, Journal of Propulsion and Power 3 (1987) 387–396. [13] W. Yang, B. Song, W. Song, Numerical study of flapping wing air vehicle based on chimera grid, in: AIAA Paper 2009-46, Presented at the 47th AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace Exposition, Orlando, Florida, Jan. 5–8, 2009. [14] H.C. Chen, E.T. Huang, Time-domain simulation of floating pier and multiple-vessel interactions by a chimera RANS method, Presented at the 7th International Symposium on fluid control, Measurement, and Visualization, Sorrento, Italy, July 25–28, 2003. [15] R.E. Gordnier, S.E. Sherer, M.R. Visbal, Computations of a maneuvering unmanned combat air vehicle using a high-order overset grid method, in: DoD High Performance Computing Modernization Program Users Group Conference, IEEE, 2007, pp. 33–40. [16] H. Pomin, S. Wagner, Aeroelastic analysis of helicopter rotor blades on deformable chimera grids, Journal of aircraft 41 (2004) 577–584.

5264

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

[17] J. Kim, D.J. Bodony, J.B. Freund, A high-order, overset-mesh algorithm for adjoint-based optimization for aeroacoustics control, in: AIAA Paper 20103818, Presented at the 16th AIAA/CEAS Aeroacoustics Conference, Stockholm, Sweden, June 7–9, 2010. [18] J. Yin, J. Delfs, Sound generation from gust-airfoil interaction using CAA-Chimera method, in: AIAA Paper 2001-2136, Presented at the 7th AIAA/CEAS Aeroacoustics Conference and Exhibit, Maastricht, Netherlands, May 28–30, 2001. [19] C.K.W. Tam, F.Q. Hu, An optimized multi-dimensional interpolation scheme for computational aeroacoustics applications using overset grids, in: AIAA Paper 2004-2812, Presented at the 10th AIAA/CEAS Aeroacoustics Conference , Manchester, Great Britain, May 10–12, 2004. [20] J. Chicheportiche, X. Gloerfelt, Direct noise computations using overset grids, in: AIAA Paper 2009-3403, Presented at the 15th AIAA/CEAS Aeroacoustics Conference, Miami, Florida, May 11–13, 2009. [21] S. Sherer, R. Gordnier, M. Visbal Gordnier, Computational study of Reynolds number and angle-of-attack effects on a 1303 UCAV configuration with a high-order overset-grid algorithm, in: Presented at the 27th AIAA Applied Aerodynamics Conference, San Antonio, Texas, June 22–25, 2009. [22] T.K. Sengupta, V.K. Suman, N. Singh, Solving Navier–Stokes equation for flow past cylinders using single-block structured and overset grids, Journal of Computational Physics 229 (2010) 178–199. [23] D.J. Bodony, G. Zagaris, A. Reichert, Q. Zhang, Provably stable overset grid methods for computational aeroacoustics, Journal of Sound and Vibration 330 (2011) 4161–4179. [24] G. Chesshire, W.D. Henshaw, Composite overlapping meshes for the solution of partial differential equations, Journal of Computational Physics 90 (1990) 1–64. [25] R. Maple, D. Belk, A new approach to domain decomposition, the Beggar code, Numerical Grid Generation in Computational Fluid Dynamics and Related Fields (1994) 305–314. [26] N.A. Petersson, Hole-cutting for three-dimensional overlapping grids, SIAM Journal on Scientific Computing 21 (1999) 646–665. [27] N.A. Petersson, An algorithm for assembling overlapping grid systems, SIAM Journal on Scientific Computing 20 (1999) 1995–2022. [28] P.G. Buning, I.T. Chiu, S. Obayashi, Y.M. Rizk, J.L. Steger, Numerical simulation of the integrated space shuttle vehicle in ascent, in: AIAA Paper 19884359, Presented at the AIAA Atmospheric Flight Mechanics Conference, Minneapolis, MN, Aug 15–17, 1988. [29] C. Kiris, D. Kwak, W. Chan, Parallel unsteady turbopump simulations for liquid rocket engines, in: Presented at the ACM/IEEE Supercomputing 2000 Conference, Dallas, Texas, November 4–10, 2000. [30] D. Brown, W. Henshaw, D. Quinlan, Overture: an object-oriented framework for solving partial differential equations, in: Y. Ishikawa, R. Oldehoeft, J. Reynders, J. Tholburn (Eds.), Scientific Computing in Object-Oriented Parallel Environments, Springer Berlin/Heidelberg, 1997, pp. 177–184. [31] W. Henshaw, Overture: an object-oriented framework for overlapping grid applications, AIAA Paper 2002-3189, Presented at the 32nd AIAA Fluid Dynamics Conference and Exhibit, St. Louis, Missouri, June 24–26, 2002. [32] Z.J. Wang, V. Parthasarathy, A fully automated Chimera methodology for multiple moving body problems, International Journal for Numerical Methods in Fluids 33 (2000) 919–938. [33] M.S. Liou, K.H. Kao, Progress in grid generation: from Chimera to DRAGON grids, Technical Report 19970031503, NASA, 1994. [34] T.L. Holst, Chimera donor cell search algorithm suitable for solving the full potential equation, Journal of Aircraft 37 (2000) 76–84. [35] J. Cai, H.M. Tsai, F. Liu, An overset grid solver for viscous computations with multigrid and parallel computing, in: AIAA Paper 2003-4232, Presented at the 16th AIAA Computational Fluid Dynamics Conference, Orlando, Florida, June 23–26, 2003. [36] B. Landmann, M. Montagnac, A highly automated parallel Chimera method for overset grids based on the implicit hole cutting technique, International Journal for Numerical Methods in Fluids (2010). [37] M. Pandya, N. Frink, R. Noack, Progress toward overset-grid moving body capability for USM3D unstructured flow solver, in: AIAA Paper 2005-5118, Presented at the 17th AIAA Computational Fluid Dynamics Conference, Toronto, Ontario, June 6–9, 2005. [38] S.E. Sherer, J.N. Scott, High-order compact finite-difference methods on general overset grids, Journal of Computational Physics 210 (2005) 459–496. [39] W.D. Henshaw, D.W. Schwendeman, An adaptive numerical scheme for high-speed reactive flow on overlapping grids, Journal of Computational Physics 191 (2003) 420–447. [40] T.L. Holst, T.H. Pulliam, Optimization of overset solution adaptive grids for hovering rotorcraft flows, in: AIAA Paper 2009-3519, Presented at the 27th AIAA Applied Aerodynamics Conference, San Antonio, Texas, June 22–25, 2009. [41] P.K. Moore, J.E. Flaherty, Adaptive local overlapping grid methods for parabolic systems in two space dimensions, Journal of Computational Physics 98 (1992) 54–63. [42] R.L. Meakin, Adaptive spatial partitioning and refinement for overset structured grids, Computer Methods in Applied Mechanics and Engineering 189 (2000) 1077–1117. [43] O. Saunier, C. Benoit, G. Jeanfaivre, A. Lerat, Third-order Cartesian overset mesh adaptation method for solving steady compressible flows, International Journal for Numerical Methods in Fluids 57 (2008) 811–838. [44] Z.J. Wang, A fully conservative interface algorithm for overlapped grids, Journal of Computational Physics 122 (1995) 96–106. [45] G. Chesshire, W.D. Henshaw, A scheme for conservative interpolation on overlapping grids, SIAM Journal on Scientific Computing 15 (1994) 819–845. [46] E. Pärt-Enander, B. Sjögreen, Conservative and non-conservative interpolation between overlapping grids for finite volume solutions of hyperbolic problems, Computers & Fluids 23 (1994) 551–574. [47] K. Mattsson, F. Ham, G. Iaccarino, Stable and accurate wave-propagation in discontinuous media, Journal of Computational Physics 227 (2008) 8753– 8767. [48] K. Mattsson, J. Nordström, High order finite difference methods for wave propagation in discontinuous media, Journal of Computational Physics 220 (2006) 249–269. [49] R.M.J. Kramer, C. Pantano, D. Pullin, Nondissipative and energy-stable high-order finite-difference interface schemes for 2-d patch-refined grids, Journal of Computational Physics 228 (2009) 5280–5297. [50] K. Mattsson, M.H. Carpenter, Stable and accurate interpolation operators for high-order multi-block finite-difference methods, SIAM Journal of Scientific Computing 32 (2010) 2298–2320. [51] J. Nordström, J. Gong, A stable hybrid method for hyperbolic problems, Journal of Computational Physics 212 (2006) 436–453. [52] J. Nordström, J. Gong, E. van der Weide, M. Svärd, A stable and conservative high order multi-block method for the compressible Navier–Stokes equations, Journal of Computational Physics 228 (2009) 9020–9035. [53] F. Olsson, N.A. Petersson, Stability of interpolation on overlapping grids, Computers & Fluids 25 (1996) 583–605. [54] G. Starius, On composite mesh difference methods for hyperbolic differential equations, Numerische Mathematik 35 (1980) 241–255. [55] J. Nordström, M. Svärd, Well-posed boundary conditions for the Navier–Stokes equations, SIAM Journal on Numerical Analysis 43 (2005) 1231–1255. [56] G.H. Golub, C.F. Van Loan, Matrix Computations, 3rd ed., Johns Hopkins University Press, 1996. [57] M.H. Carpenter, J. Nordström, D. Gottlieb, A stable and conservative interface treatment of arbitrary spatial accuracy, Journal of Computational Physics 148 (1999) 341–365. [58] 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, Journal of Computational Physics 111 (1994) 220–236. [59] K. Mattsson, J. Nordström, Summation by parts operators for finite difference approximations of second derivatives, Journal of Computational Physics 199 (2004) 503–540. [60] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time Dependent Problems and Difference Methods, Wiley, 1995. [61] B. Strand, Summation by parts for finite difference approximations for d/dx, Journal of Computational Physics 110 (1994) 47–67. [62] H.-O. Kreiss, G. Scherer, Finite element and finite difference methods for hyperbolic partial differential equations, in: C. De Boor (Ed.), Mathematical Aspects of Finite Elements in Partial Differential Equations, Academic Press, 1974, pp. 195–212.

A. Reichert et al. / Journal of Computational Physics 231 (2012) 5243–5265

5265

[63] A.E. Chertock, Strict stability of high-order compact implicit finite-difference schemes: the role of boundary conditions for hyperbolic PDEs, Ph.D. Thesis, Tel-Aviv University, 1998. [64] S.S. Abarbanel, A.E. Chertock, Strict stability of high-order compact implicit finite-difference schemes: The role of boundary conditions for hyperbolic PDEs, I, Journal of Computational Physics 160 (2000) 42–66. [65] S.S. Abarbanel, A.E. Chertock, A. Yefet, Strict stability of high-order compact implicit finite-difference schemes: The role of boundary conditions for hyperbolic PDEs, II, Journal of Computational Physics 160 (2000) 67–87. [66] M. Svärd, J. Nordström, A stable high-order finite difference scheme for the compressible Navier–Stokes equations: no-slip wall boundary conditions, Journal of Computational Physics 227 (2008) 4805–4824. [67] M. Svärd, M.H. Carpenter, J. Nordström, A stable high-order finite difference scheme for the compressible Navier–Stokes equations, far-field boundary conditions, Journal of Computational Physics 225 (2007) 1020–1038. [68] P. DuChateau, D. Zachman, Applied Partial Differential Equations, Dover, New York, USA, 2002. [69] M. Svärd, On coordinate transformations for summation-by-parts operators, J. Sci. Comput. 20 (2004) 29–42. [70] J. Nordström, Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation, J. Sci. Comput. 29 (2006) 375–404. [71] A. Reichert, Stable numerical methods for hyperbolic partial differential equations using overlapping domain decomposition, Ph.D. Thesis, University of Illinois at Urbana-Champaign, 2011. [72] K. Mattsson, M. Svärd, M. Carpenter, J. Nordström, High-order accurate computations for unsteady aerodynamics, Computers & Fluids 36 (2007) 636– 649. [73] S.K. Lele, Compact finite difference schemes with spectral-like resolution, Journal of Computational Physics 103 (1992) 16–42. [74] R. Vichnevetsky, J.B. Bowles, Fourier Analysis of Numerical Approximations of Hyperbolic Equations, Society for Industrial Mathematics, 1982.