Summation-by-Parts Operators for Non-Simply Connected Domains

(1)

Vol. 40, No. 3, pp. A1250–A1273

SUMMATION-BY-PARTS OPERATORS FOR NON-SIMPLY

CONNECTED DOMAINS∗

SAMIRA NIKKAR† AND JAN NORDSTR ¨OM†

Abstract. We construct fully discrete stable and accurate numerical schemes for solving partial differential equations posed on non-simply connected spatial domains. The schemes are constructed using summation-by-parts operators in combination with a weak imposition of initial and bound-ary conditions using the simultaneous approximation term technique. In the theoretical part, we consider the two-dimensional constant coefficient advection equation posed on a rectangular spatial domain with a hole. We construct the new scheme and study well-posedness and stability. Once the theoretical development is done, the technique is extended to more complex non-simply connected geometries. Numerical experiments corroborate the theoretical results and show the applicability of the new approach and its advantages over the standard multiblock technique. Finally, an application using the linearized Euler equations for sound propagation is presented.

Key words. initial boundary value problems, stability, well-posedness, boundary conditions, non-simply connected domains, complex geometries

AMS subject classification. 65M99 DOI. 10.1137/18M1163671

1. Introduction. High order summation-by-parts (SBP) operators, together with a weak imposition of initial and well-posed boundary conditions using the si-multaneous approximation term (SAT) technique, provide provably fully discrete un-conditionally stable schemes for steady or time-dependent spatial domains [1, 2, 3, 5]. These schemes have so far been mostly developed for spatial domains consisting of sim-ply connected regions. To handle more complicated geometries, hybrid formulations utilizing finite volume and finite difference methods [6, 7, 8, 9] have been proposed. Other alternatives within the finite difference community for complex geometries in-clude finite difference schemes using over-set mesh discretizations [10, 11, 12, 13], multiblock techniques [14, 15, 16, 17], as well as SBP extensions to unstructured grids [18, 19, 20].

In this article, we extend the SBP-SAT technique to handle partial differential equations posed on non-simply connected multidimensional geometries. As a proto-type problem, we consider rectangles with rectangular holes inside, where the holes are not part of the computational domain. Hybrid methods, discontinuous Galerkin methods, and other types of schemes of SBP-SAT form are easily applicable to non-simply connected domains. In these methods, one divides the spatial domain into patches and glues the numerical solutions together via interface couplings. In this paper, we proceed the other way around and discretize the geometry patch surround-ing the hole directly, without subdividsurround-ing. The novelty of the new technique is that it reduces the number of block-to-block connections and in that way improves the accuracy and lowers the cost.

The rest of the article proceeds as follows. In section 2, we study the two-dimensional constant coefficient advection equation posed on a rectangular geometry

∗_{Submitted to the journal’s Methods and Algorithms for Scientific Computing section January 3,}

2018; accepted for publication (in revised form) January 23, 2018; published electronically May 1, 2018.

http://www.siam.org/journals/sisc/40-3/M116367.html

†_{Department of Mathematics, Computational Mathematics, Link¨}_{oping University, SE-581 83}

Link¨oping, Sweden (samira.nikkar@gmail.com, jan.nordstrom@liu.se).

A1250

(2)

Fig. 1. (a) The model geometry, Ω, its boundaries denoted by A, B, C, D, a, b, c, d, and the corresponding outward pointing normal vectors from each boundary. (b) A schematic of the

different parts of δΩ influenced by the nonzero terms used in the construction of P_C,C,D,D−1 in (17).

with a hole. In section 3, the discrete problem and a new combination of SBP op-erators are presented. Stability of the new scheme is investigated in section 4. We extend the new approach to more complex geometries in sections 5 and 6. Numerical calculations are shown in section 7, where we measure the accuracy and efficiency of our scheme and compare it to the standard SBP-SAT multiblock technique. Finally, we summarize and draw conclusions in section 8.

2. Well-posedness. To develop the theory, we consider the constant coefficient scalar advection equation in two space dimensions,

(1) ut+ αux+ βuy = 0, (x, y) ∈ Ω, t ∈ [0, T ],

where the subscripts t, x, and y denote partial derivatives. The computational domain Ω is depicted in Figure 1(a). The spatial region H is not a part of the computational domain; it forms a hole in Ω.

The energy method (multiplying (1) with the solution and integrating over Ω) together with the use of the Gauss–Green theorem gives

(2) ||u||2

t = −

I

δΩ

u2(α, β) · n ds,

where δΩ = {A∪B ∪C ∪D ∪a∪b∪c∪d} is the boundary of Ω. Moreover, n = (n1, n2)

is the outward pointing normal vector from Ω and ds is an infinitesimal element along δΩ. The norm is defined as ||u||2₌RR

Ωu

2 _{dx dy.}

In order to bound the energy rate of the solution in (2), we specify

(3) u = g if (α, β) · n < 0,

where (α, β) · n = n1α + n2β. Assuming, for example, that α, β > 0, (3) leads to

(4) us= gsand uw= gw,

(3)

where s ∈ {C, c} and w ∈ {D, d}; see Figure 1(a). We insert (4) into (2), integrate in time and consider an initial condition u = f . The continuous energy estimate becomes (5) ||u(T )||2= ||f ||2 + β X s∈{C,c} Z T 0 Z s g2_s dx dt + αX w∈{D,d} Z T 0 Z w g2_w dy dt + BT.

In (5), P denotes summation and BT is the negative contribution from the outflow boundaries {A, B, a, b} as (6) BT = −β X n∈{A,a} Z T 0 Z n u2 dx dt − α X e∈{B,b} Z T 0 Z e u2dy dt.

We summarize the result in the following.

Proposition 1. The continuous problem (1) for α, β > 0 augmented with bound-ary conditions (4) is strongly well-posed and has the bound (5).

For more details on well-posed and strongly well-posed problems, see [21]. 3. Summation-by-part operators. The domain Φ = {Ω ∪ H} is a rectangle and we discretize it using N and M grid points in the x and y directions, respectively. In time we use L time levels. The boundaries of H coincide with the coordinate lines after the discretization of Φ. We allocate a column vector of size LM N to the grid as

(7) v =         v1 .. . [vk] .. . vL         , [vk] =         v1 .. . [vi] .. . vN         k , [vi]_k=         v1 .. . vj .. . vM         ki in which

(8) vkij≈ u(tk, xi, yj) for (xi, yj) ∈ Ω, vkij: not defined for (xi, yj) ∈ H.

Inside the hole vkij has neither any relation to the solution u nor any contribution

to the calculations (we will show the latter below). We have related v to these grid points for the convenience of using tensor products in the formulations below.

The first derivative uy at xi for all i ∈ {1, . . . , Nl} ∪ {N − Nr+ 1, . . . , N } (see

Figure 2) is approximated by DM

y u. In Figure 2, and also in the remainder of this

article, Nl and Nr are the number of grid points, in the x direction, on the

left-and right-hleft-and sides of the hole, respectively. Additionally, Ma and Mb denote the

number of grid points, in the y direction, above and below the hole, respectively. DMy

is a so-called SBP operator of the form

(9) DyM = (PyM)−1QMy ,

and u = [u1, . . . , uM]T is a smooth function injected in each grid point in the y

direction. The superscript M denotes the size of the matrix and the subscript y denotes the direction along which the operator is acting. Moreover, PM

y is a symmetric

positive definite matrix, and QM

y is an almost skew-symmetric matrix that satisfies

(10) QM_y + (QM_y )T= E₁M−EM

0 = B

M _{= diag(−1, 0, . . . , 0, 1).}

(4)

In (10), EM

0 = diag(1, 0, . . . , 0) and E1M= diag(0, . . . , 0, 1). The first derivative in the x

direction, DN

x = (PxN)−1QNx, at yj for all j ∈ {1, . . . , Mb} ∪ {M − Ma+ 1, . . . , M } (see

Figure 2), and the first derivative in time, DL

t = (PtL)−1QLt, are defined in analogous

ways.

The first derivative uy at xi for all i ∈ {Nl+ 1, . . . , N − Nr} is approximated by

˜ Dyu, where ˜Dy is (11) D˜y=   DMb y 0M −(Ma+Mb) DMa y  . In (11), DMb y = (P Mb y )−1Q Mb y and D Ma y = (P Ma y )−1Q Ma

y are the same type of SBP

operators as in (9) but smaller in size (Mb and Ma, respectively). The notation

0 denotes a zero matrix and the superscript denotes its size (this notation is used throughout the rest of the paper). The zero matrix corresponds to the grid points inside the hole. The derivative in the x direction at yjfor all j ∈ {Mb+1, . . . , M −Ma}

is constructed in the same way, as

˜ Dx=   DNl x 0N −(Nl+Nr) DNr x  . (12)

A finite difference approximation including the time discretization [22, 23], on SBP form, is constructed by extending the one-dimensional SBP operators in a tensor product fashion as Dt= DLt ⊗ [I_xΩ⊗ IΩ y] + [I H x ⊗ I Ω y] + [I Ω x ⊗ I H y ] , (13) Dx= It⊗ [DN_x ⊗ IΩ y] + [ ˜Dx⊗ IyH] , Dy = It⊗ [I_xΩ⊗ DM y ] + [I H x ⊗ ˜Dy] ,

where Itis the identity matrix in time and has the size L. In (13), ⊗ represents the

Kronecker product which is defined as

(14) A ⊗ B =      a11B a12B . . . a1nB a21B a22B . . . a2nB .. . ... ... am1B am2B . . . amnB     

for arbitrary matrices A and B. For an m × n matrix A and a k × l matrix B, the size of A ⊗ B is (mk) × (nl). More details on Kronecker products and their properties can be found in [24, 25].

Moreover, IH

y and IxH are diagonal matrices of size M and N respectively given

by (15) IH y =   0Mb IM −(Ma+Mb) 0Ma  and IxH =   0Nl IN −(Nl+Nr) 0Nr  ,

where I (with a slight abuse of notation) denotes the identity matrix and the super-scripts the size of the matrices. (In other words, IH

x and IyH have diagonal elements

(5)

x i=1 i=Nl i=N-‐Nr+1 i=N

𝑦 j=1 j=Mb j=M-‐Ma+1 j=M 𝐷!! 𝐷!! 𝐷!! 𝐷!! 𝐷!!! 𝐷!!! 𝐷!!! 𝐷!!!

Fig. 2. A schematic of the regions where difference operators in the x and y directions are defined. x y 𝐼_!! 𝐼!! 𝐼!! 𝐼!!⊗ 𝐼!! 𝐼!!⊗ 𝐼!! _𝐼_!!_{⊗ 𝐼} !! 𝐼!!⊗ 𝐼!! 𝐼!!⊗ 𝐼!! 𝐼!!⊗ 𝐼!! 𝐼!!⊗ 𝐼!! 𝐼!!⊗ 𝐼!! 𝐼!! 𝐼_!! 𝐼!!

Fig. 3. A schematic of the geometry

in-fluenced with different combinations of IΩ

x,yand IH x,y. x y 𝐼!!⨂𝐷!! 𝐼!!⊗ 𝐷!! 𝐼!! 𝐼!! 𝐼!!⨂𝐷!! 𝐼_!!_{⊗ 𝐷}!_! 𝐷_!! 𝐷!_!

Fig. 4. A schematic that shows where the

matrices, DN

x, D˜x, IyΩ, and IyH, are defined;

the solid lines correspond to the nonzero con-tributions. x y 𝐷!!⊗ 𝐼!! 𝐼!! 𝐼!! 𝐷_!!_{⊗ 𝐼} !! 𝐷!_! _{⨂ 𝐼} ! ! 𝐷!! 𝐷!_! 𝐷!! ⨂ 𝐼!!

matrices, DM

y , D˜y, IxΩ, and IxH, are defined;

the solid lines correspond to the nonzero con-tributions.

equal to one corresponding to the points inside the hole, respectively, and zero diag-onal elements otherwise.) Moreover, IxΩ= IN − IxH and IyΩ= IM− IyH, where IM,N

are identity matrices of size M and N , respectively.

A schematic of different regions of Ω influenced by the different identity matrices IΩ

x,y and Ix,yH is shown in Figure 3. Additionally, schematics of the different regions

where the matrices Dxand Dyare active, as well as the underlying matrices involved

in their construction, are shown in Figures 4 and 5.

Lemma 2. The spatial difference operators in (13) commute. To see the proof, see [3].

(6)

4. Stability. The fully discrete SBP-SAT approximation of (1), including (4), can be written as Dtv+αDxv+βDyv = Pi−1σi(vi−f )+ X s={C, c} P_s−1σs(vs−gs)+ X w={D, d} P_w−1σw(vw−gw), (16)

where σi,s,ware penalty coefficients for the weak initial and boundary conditions. gs,w

are zero vectors of the same size as v, except at the positions corresponding to the inflow boundaries where the zeros are replaced with the boundary data. Moreover f is a zero vector, of the same size as v, except at the positions corresponding to t = 0 where the initial data (compatible with the boundary conditions) is injected. The subscripts i, s, and w on the solution restrict the solution to the initial time, and the s and w boundary locations. Additionally,

P_i−1= (P_tL)−1E₀L⊗ [IxΩ⊗ I Ω y] + [I H x ⊗ I Ω y] + [I Ω x ⊗ I H y ] , (17) P_C−1= It⊗ [I_xΩ⊗ (PyM)−1E M 0 ] + [I H x ⊗ ( ˜Py)−1E˜0b] , P_c−1= It⊗ IxH⊗ ( ˜P ) −1 y E˜0a, P_D−1= It⊗ [(P_xN)−1E₀N ⊗ IΩ y] + [( ˜Px)−1E˜0l⊗ I H y ] , P_d−1= It⊗ ( ˜Px)−1E˜0r⊗ I H y , where (18) ( ˜Py)−1=   (PMb y )−1 0M −(Ma+Mb) (PMa y ) −1  , ( ˜Px)−1=   (PNl x )−1 0N −(Nl+Nr) (PNr x ) −1  , ˜ E0b=   EMb 0 0M −(Ma+Mb) 0Ma  , E˜0a=   0Mb 0M −(Ma+Mb) EMa 0  , ˜ E0r=   0Mb 0M −(Ma+Mb) ENr 0  , E˜0l=   ENl 0 0N −(Nl+Nr) 0Nr  .

Note that we have again slightly abused notation in (18) by applying the inverse sign on ˜Pyand ˜Px. In Figure 1(b), we graphically show the different parts of δΩ (the

inflow parts) that are influenced by nonzero contributions from the penalty terms in (17).

Next, we apply the discrete energy method, by multiplying (16) from the left with vT_{P , where} (19) P = PtL⊗ [PxN ⊗ IyΩ] + [ ˜Px⊗ IyH] | {z } :=Px [IxΩ⊗ PyM] + [IxH⊗ ˜Py] | {z } :=Py , (20) ˜ Py=   PMb y 0M −(Ma+Mb) PMa y   and P˜x=   PNl x 0N −(Nl+Nr) PNr x  . Schematics of the different regions in the geometry influenced by Px and Py, as well

as the underlying matrices involved in their construction, are shown in Figures 6 and 7. By using the properties of the Kronecker product [24, 25] one can rewrite (19) as

(7)

x y 𝐼!!⨂𝑃!! 𝐼!!⊗ 𝑃!! 𝐼!! 𝐼!! 𝐼!!⨂𝑃!! 𝐼!!⊗ 𝑃!! 𝑃!! 𝑃!!

matrices, PM

y , P˜y, IxΩ, and IxH, are defined;

the solid lines correspond to the nonzero con-tributions. x y 𝑃!!⊗ 𝐼!! 𝑃!! ⨂ 𝐼!! 𝐼!! 𝐼!! 𝑃!!⊗ 𝐼!! 𝑃!! ⨂ 𝐼!! 𝑃!! 𝑃!!

matrices, PN

x , P˜x, IΩy, and IyH, are defined;

the solid lines correspond to the nonzero con-tributions. x y 𝑃! 𝑃! 𝑃! 𝑃! 𝑃! 𝑃! 𝑃! 𝑃!

Fig. 8. A schematic of the geometry and its relation to different parts of the matrix P .

(21) P = P_tL⊗   I Ω xP N x ⊗ I Ω yP M y | {z } :=P1 + I_xHP_xN ⊗ IΩ yP˜y | {z } :=P2 + I_xΩP˜x⊗ IyHP M y | {z } :=P3   .

In Figure 8, the nonzero contributions from P1,2,3to the regions of Ω are shown.

We need the following lemma.

Lemma 3. The matrix P in (21) is nonsingular. Proof. Consider the matrix

(22) S = (PtL)−1⊗ S1+ S2+ S3

(8)

in which S1= IxΩ(P N x ) −1_{⊗ I}Ω y(P M y ) −1_, (23) S2= IxH(P N x )−1⊗ I Ω y( ˜Py)−1, S3= IxΩ( ˜Px)−1⊗ IyH(P M y ) −1_.

The properties of the Kronecker product give

(24) SP = It⊗ [(S1+ S2+ S3)(P1+ P2+ P3)].

Now, we substitute S1,2,3 from (23) and P1,2,3 from (21) in (24) and find

(S1+ S2+ S3)(P1+ P2+ P3) = S1P1+ S2P2+ S3P3 = [I_xΩ⊗ IyΩ] + [I H x ⊗ I Ω y] + [I Ω x ⊗ I H y ] := I Ω

by the fact that SiPj = 0 if i 6= j. In the same way, we find that P S = It⊗ IΩ.

Applying the discrete energy method to (16) and considering zero data gives

vTP Dtv+αvTP Dxv+βvTP Dyv = vTP Pi−1σivi+vTP PC−1σCvC+vTP Pc−1σcvc

(25)

+ vTP P_D−1σDvD+vTP P_d−1σdvd.

Next, we evaluate the matrix products in (25) as follows: P Dt= QLt ⊗ (P1+ P2+ P3), P Dx= PtL⊗ [I_xΩQN_x ⊗ IyΩP M y ] + [I H x Q N x ⊗ I Ω yP˜y] + [IxΩQ˜x⊗ IyHP M y ] , P Dy = PtL⊗ [I_xΩP_xN ⊗ IyΩQ M y ] + [I Ω xP˜x⊗ IyHQ M y ] + [I H x P N x ⊗ I Ω yQ˜y] , (26) where (27) Q˜x=   QNl x 0N −(Nl+Nr) QNr x   and Q˜y=   QMb y 0M −(Ma+Mb) QMa y  . Further, (28) P P_i−1= E0L⊗ (P1+P2+P3) := Hi, P P_C−1= PL t ⊗ [IΩ xPxN ⊗ IyΩEM0 ] + [IxHPxN ⊗ IyΩE˜0b] := HC, P Pc−1= PtL⊗ IxHPxN ⊗ IyΩE˜0a:= Hc, P P_D−1= PtL⊗ [IxΩE0N ⊗ IyΩPyM] + [IxΩE˜0l⊗ I H y PyM] := HD, P P_d−1= PL t ⊗ IxΩE˜0r⊗ I H y PyM := Hd.

The details of the computations in (26) and (28) are given in Appendix A. By substituting (26) and (28) into (25) and adding the transpose, we obtain

vTHfv − vT(1 + 2σi)Hiv = (β + 2σc)vTHcv + (β + 2σC)vTHCv

(29)

+ (α + 2σd)vTHdv + (α + 2σD)vTHDv

+ CT,

(9)

SBP21, x derivative SBP21, y derivative SBP42, x derivative SBP42, y derivative

Fig. 9. A closer look at CT , near one of the corners, for SBP 21 and SBP 42 schemes. Table 1

The size of each block in CT around one corner for schemes of different order of accuracy.

SBP 21 42 63

Size of CT 2 × 1 4 × 4 6 × 6

where Hf = E1L⊗ (P1+ P2+ P3) and CT stands for corner terms. Details of the

derivation of (29) are given in Appendix B.

In (29), we choose σi= −1, σc,C = −β, and σd,D= −α which, in the absence of

CT and the presence of data, lead to

(30) ||v||2Hf=||f || 2 Hi+ β X s∈{C,c} ||gs||2Hs+α X w∈{D,d} ||gw||2Hw+DI,

where DI includes the extra damping effects from the initial and inflow boundary procedures. However, in the presence of CT , on the right-hand side of (30) we have an indefinite term that involves the solution on a few grid points enclosed in two small blocks around each corner. These blocks are results of using central differences in the ξ and η directions around the corners while having different norms in the perpendicular directions. Hence, the skew-symmetric property of the difference operators cannot be preserved and the interior grid contribution is not removed in CT . In Figure 9 we show schematically where CT is located when using the SBP21 and SBP42 operators. The size of each block depends on the order of the difference operators and is independent of the number of grid points; see Table 1.

In order to investigate the stability, we consider the semidiscrete version of (1) in SBP-SAT form, written as

(31) vt+ Avx= 0,

(10)

0 20 40 60 80 100 120 Real -200 -150 -100 -50 0 50 100 150 200 Imaginary SBP21, N=M=91

Fig. 10. The discrete spectrum for the sec-ond order case (SBP 21), N = M = 91.

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Real -150 -100 -50 0 50 100 150 Imaginary SBP21, N=M=91

Fig. 11. A blow-up of the spectrum near the imaginary axis, SBP 21, N = M = 91.

0 10 20 30 40 50 60 Real -250 -200 -150 -100 -50 0 50 100 150 200 250 Imaginary SBP42, N=M=91

Fig. 12. The discrete spectrum for the

third order case (SBP 42), N = M = 91.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Real -250 -200 -150 -100 -50 0 50 100 150 200 250 Imaginary SBP42, N=M=91

and investigate the spectrum of A. In (31), A is given by

(32) A = αDx+ βDy− X s={C,c} P_s−1σs− X w={D,d} P_w−1σw,

and σs,wand PC,c,D,d−1 are given in (16) (discard Itin (17)). We assume α = 1, β = 1

and use the following stability conditions for the penalty parameters: (33) 1 + 2σi≤ 0, β + 2σC,c≤ 0, α + 2σD,d≤ 0.

We choose σi= −1, σC= σc = −β, and σD= σd= −α.

The eigenvalue distribution of A with SBP operators of different orders on a grid of size 91 × 91 are shown in Figures 10–17. The minimum real part of the spectrum for a sequence of mesh refinements is given in Figure 18. As Figures 10–18 show, the eigenvalues of A for all orders of accuracy have the correct sign with a minimum real part clearly positive. In combination with SBP-SAT in time, this implies stability for these particular meshes.

5. Extension to geometries with multiple holes. One can readily extend the techniques presented in sections 3 and 4 to construct SBP operators and stable schemes for geometries with multiple holes; see Figure 19. In this example, we parti-tion Ω along the x and y axes, as seen in Figure 20. The partiparti-tioning along the x axis

(11)

0 10 20 30 40 50 60 Real -300 -200 -100 0 100 200 300 Imaginary SBP63, N=M=91

Fig. 14. The discrete spectrum for the

fourth order case (SBP 63), N = M = 91.

0 0.05 0.1 0.15 0.2 Real -300 -200 -100 0 100 200 300 Imaginary SBP63, N=M=91

0 10 20 30 40 50 60 Real -400 -300 -200 -100 0 100 200 300 400 Imaginary SBP84, N=M=91

Fig. 16. The discrete spectrum for the fifth order case (SBP 84), N = M = 91. 0 0.02 0.04 0.06 0.08 0.1 Real -400 -300 -200 -100 0 100 200 300 400 Imaginary SBP84, N=M=91

60 80 100 120 0 0.02 0.04 0.06 0.08 0.1

Number of grid points in each space dimension

min(real( λ )) SBP84 SBP42 SBP63 SBP21

Fig. 18. The minimum real part of the spectrum for a sequence of mesh refinements for different orders of accuracy.

is done such that in each partition, only one difference operator approximates the y derivative. The partitioning along the y axis is analogous.

The difference operators in the y and x directions are defined using tensor products in the following way:

(34) Dy = 5 X i=1 (INi x ⊗ D i y), Dx= 5 X j=1 (D_xj⊗ IMj y ).

(12)

𝑦

x

𝐻

_!

𝐻

_!

Ω

Fig. 19. A geometry with two holes. x 𝑦 𝑦 x 𝐷!!,!,! 𝐷!! 𝐷!! 𝐷!!,! 𝐷!! 𝐷!! 𝐷!!

Fig. 20. A schematic of the partitioning along the x and y axes.

In (34), the matrices INi

x and I

Mj

y are defined to single out the different segments

on the x and y axes over which one partition is defined. As an example, IN1

x is an

N × N matrix (N =P5

i=1Ni) which has elements equal to one on the main diagonal

corresponding to the N1 grid points shown in Figure 21. Other matrices are defined

similarly. Moreover, in (34), we have used

D1_y= (P_yM)−1QM_y , (35) D2_y=   (PM1+M2+M3 y )−1QMy 1+M2+M3 0M4 (PM5 y )−1QMy5  ,

(13)

A1262 SAMIRA NIKKAR AND JAN NORDSTR ¨OM

𝐼

_!!!

𝐼

_!!!

𝐼

_!!!

𝐼

_!!!

𝐼

_!!!

x

𝑦

𝐼

_!!!

𝐼

_!!!

𝐼

_!!!

𝐼

_!!!

𝐼

_!!!

Fig. 21. A schematic of the different segments along the x and y axes influenced by

IN1,N2,N3,N4,N5

x and IyM1,M2,M3,M4,M5, respectively. The superscripts denote the number of grid

points in each segment; the solid lines correspond to the nonzero elements; end points are boldface if they belong to that term.

D3_y=       (PM1 y )−1QMy1 0M2 (PM3 y )−1QMy3 0M4 (PM5 y )−1QMy 5       , D4_y=   (PM1 y )−1QMy1 0M2 (PM3+M4+M5 y )−1QMy3+M4+M5  , D5y= Dy1= (PyM)−1QMy . Similarly we have D1_x= (P_xN)−1QN_x, (36) D2_x=   (PN1+N2 x )−1Q N1+N2 x 0N3+N4 (PN5 x )−1QNx5  , D3_x= D1_x = (P_xN)−1QN_x, D4_x=   (PN1 x )−1QNx1 0N2+N3 (PN4+N5 x )−1QNx4+N5  , D5_x= D1_x = (P_xN)−1QN_x.

(14)

𝐻

_!

Ω

x

y

ξ

𝜂

Ω

!

𝐻!

!

Fig. 22. A more general non-simply connected geometry in Cartesian coordinates transformed into body-fitted coordinates.

The SBP-SAT scheme is similar but more technically involved than the ones in section 3.

6. Extension to curved geometries with holes. We extend the well-posedness analysis, SBP operators constructions, and the stability procedure pre-sented in sections 2, 3, and 4 to a general geometry with a hole. As seen in Figure 22, the curved computational domain is mapped into a rectangular geometry with a hole, namely, ˜Ω, by using a Cartesian to the body-fitted coordinate transformation (x, y) *) (ξ, η). The Jacobian matrix of the transformation is denoted by [J ], where

(37) [J ] =

xξ yξ

xη yη

and (∂/∂x, ∂/∂y)T = [J ](∂/∂ξ, ∂/∂η)T. Equation (1) is then transformed to the body-fitted coordinate system, which results in a variable coefficient advection equa-tion,

(38) J ut+ ˜α(ξ, η)uξ+ ˜β(ξ, η)uη= 0, (ξ, η) ∈ ˜Ω, t ∈ [0, T ].

In (38), J = xξyη − xηyξ is the determinant of [J], ˜α(ξ, η) = αJ ξx + βJ ξy, and

˜

β(ξ, η) = αJ ηx+ βJ ηy, where ˜α, ˜β are variables and α, β are constants.

The energy method applied to (38) together with the use of the Green–Gauss theorem in ˜Ω gives (39) d dt(||u|| 2 J) = − I δ ˜Ω u2( ˜α, ˜β) · n ds,

where δ ˜Ω = {A∪B ∪C ∪D ∪a∪b∪c∪d} is the boundary of ˜Ω. Moreover, n = (n1, n2)

is the outward pointing normal vector from ˜Ω and ds is an infinitesimal element along δ ˜Ω. The norm is defined as ||u||2

J=

RR

˜

ΩJ u

2 _{dξ dη.}

In order to bound the energy rate of the solution in (39), we specify

(40) u = g if ( ˜α, ˜β) · n < 0,

where ( ˜α, ˜β) · n = n1α(ξ, η) + n˜ 2β(ξ, η). For clarity and brevity, we only consider˜

the s and w boundaries (see Figure 1(a)), where s ∈ {C, c} and w ∈ {D, d}. The boundary conditions are chosen to be

(15)

us= gs if β˜s> 0,

(41)

uw= gw if α˜w> 0.

We insert (41) into (39), integrate in time, and consider an initial condition u = f . The continuous energy estimate becomes

||u(T )||2J= ||f || 2 J+ ˜ βs+ | ˜βs| 2 X s∈{C,c} Z T 0 Z s g2_s dξ dt (42) +α˜w+ | ˜αw| 2 X w∈{D,d} Z T 0 Z w g2_wdη dt + ˜ βs− | ˜βs| 2 X s∈{C,c} Z T 0 Z s u2s dξ dt +α˜w− | ˜αw| 2 X w∈{D,d} Z T 0 Z w u2_w dη dt + BT.

In (5), BT is the contribution from the other boundaries {A, B, a, b} which are treated similarly to the s and w boundaries.

We discterize the transformed geometry using the standard technique described in section 3. However, the governing equation needs to be split [3, 4] prior to the discretization. The split version of (38) is

(43) J ut+ ( ˜αu)ξ− ˜αξu + ( ˜βu)η− ˜βηu = 0, (ξ, η) ∈ ˜Ω, t ∈ [0, T ].

The fully discrete SBP-SAT approximation of (43), including (41), is DtJ v + DξAv − Aξv + DηBv − Bηv (44) = P_i−1σi(vi− f ) + X s={C, c} P_s−1Σs(vs− gs) + X w={D, d} P_w−1Σw(vw− gw),

where σi is the initial penalty coefficient and Σs,w are penalty matrices for the weak

boundary conditions. Moreover, gs,ware vectors of the same size as v, with boundary

data at appropriate positions and zeros elsewhere. Additionally, f is a zero vector, of the same size as v, except at the positions corresponding to t = 0 where the initial data is injected. The subscripts i, s, and w on the solution restrict the solution to the initial time and the s and w boundary locations. The matrices P_i,s,w−1 are defined similarly to their counterparts in (17). Finally, J , A, and B are diagonal matrices approximating J , ˜α, and ˜β, defined by

(45) Aξ = Dξ(DηηA − DηξB), Bη= Dη(−DξηA + DξξB).

In (45), ξ and η are diagonal matrices having the ξ and η coordinates at the relevant locations on the diagonal.

The discrete energy method with zero data result in

vTP J DtJ v+vTP DξAv+vTP DyBv+vTP (Aξ+ Bη)v = vTP Pi−1σivi

+ vTP P_C−1ΣCvC+vTP Pc−1Σcvc+vTP PD−1ΣDvD+vTP Pd−1Σdvd.

(46)

(16)

Following a similar procedure as in the derivation of (26)–(28) one obtains vTHfv + vTP (Aξ+ Bη)v = vT(B + 2Σc)Hcv + vT(B + 2ΣC)HCv

(47)

+ vT(A + 2Σd)Hdv + vT(A + 2ΣD)HDv

+ vT(1 + 2σi)Hiv + CT,

where Hf = J (E1L⊗ (P1+ P2+ P3)) and CT stands for corner terms.

As a result of Lemma 2, Aξ+ Bη = 0 (the numerical geometric conservation law

holds) and hence the last term on the left-hand side of (47) vanishes.

To get an estimate as in (30), we choose σi = −1, Σc,C = −(|B| + B/2), and

Σd,D = −(|A| + A)/2. These choices together with nonzero data and disregarding

CT lead to (48) ||v||2 Hf = ||f || 2 Hi+ B+|B| 2 P s||gs||2Hs+ A+|A| 2 P w||gw||2Hw+ DI, where DI is the extra dissipation.

Remark 4. The extension to general curvilinear geometries with multiple holes of various kinds can be handled with a generalization of the technique discussed in section 5.

7. Numerical experiments. We consider the two-dimensional constant coeffi-cient symmetrized Euler equations [26]

(49) Ut+ ˆAUx+ ˆBUy= 0, (x, y) ∈ Ω, t ∈ [0, T ],

where U = [¯cρ/√γ ¯ρ, u, v, θ/(¯cpγ(γ − 1))]T_{. In (49), ρ, u, v, θ, and γ are the}

density, the x and y velocity components, the temperature, and the ratio of specific heats, respectively. An equation of state of the form γp = ¯ρθ + ρ¯θ, where p is the pressure, closes the system (49). Moreover, the bar sign denotes the state around which we have linearized. The matrices in (49) are

(50) A =ˆ       ¯ u c/¯√γ 0 0 ¯ c/√γ u¯ 0 qγ−1_γ c¯ 0 0 u¯ 0 0 qγ−1_γ c¯ 0 u¯       , ˆB =       ¯ v 0 ¯c/√γ 0 0 ¯v 0 0 ¯ c/√γ 0 v¯ q γ−1 γ ¯c 0 0 qγ−1_γ ¯c v¯       .

In the remainder of this paper, we use γ = 1.4, ¯c = 2, and ¯ρ = 1.

7.1. Accuracy. To verify the order of accuracy of the scheme we use the domain in Figure 1(a), where Φ = [0, 1] × [0, 1] and H = [1/3, 2/3] × [1/3, 2/3]. We prescribe the mean velocity field to be (¯u, ¯v) = (1, 1) and use the manufactured solution (51) U∞= [sin(x − t), cos(x − t), sin(y − t), cos(y − t)]T

for the forcing function and initial and boundary data in (49). Moreover, characteristic boundary conditions [3] are used.

We examine the scheme for SBP operators of order 2s in the interior and s close to the boundaries in space for s ∈ {1, 2, 3, 4}. The fifth order accurate SBP operator, with L = 81 and T = 0.1, is used in time. The rates of convergence are shown in Figure 23. According to [27, 28, 29, 5, 30], for a scheme with first derivative SBP operators which are 2s order accurate in the interior and s order accurate close to the boundaries (where diagonal norms are used) should yield s + 1 order of accuracy globally. The results in Figure 23 converge at these rates.

(17)

102

Number of grid points in each space dimension 10-10 10-8 10-6 10-4 The error SBP21 s=2 SBP42 s=3 SBP63 s=4 SBP84 s=5

Fig. 23. Mesh refinement, the errors, and the convergence rates.

x y Interface Interface Interface Interface Interface Interface Interface Interface

Fig. 24. A schematic of the

multi-domains and interfaces.

7.2. Comparison with the multiblock technique. In the standard SBP-SAT multiblock approach, the domain is divided into subdomains and interface penal-ties are used to couple the blocks. A schematic of such multiblock division including the interfaces (dashed lines) is shown in Figure 24.

One drawback with the standard multiblock schemes is that close to interfaces, the accuracy of the approximation is reduced from 2s to s. In order to compare this approach with our new one, we use the domain in Figure 1(a). To compare these schemes we consider the manufactured solution

(52) U∞= [0, 0, 0, e−10((x−t)

2_+(y−t)2₎

]T.

The resulting solution is a pressure pulse that starts from (x, y) = (0, 0) and travels to (x, y) = (1, 1) as time passes.

For low order operators and coarse meshes the influence of the interfaces is visible on the solution, as seen in Figures 25–28, where the pressure distribution for different times is shown. We have used SBP 21 and a coarse grid of size 19 × 19 in space together with a fifth order SBP operator (with sufficiently small time steps) in time.

The interfaces are clearly visible in the multiblock approach, while in our approach they are not. For finer grids and higher order operators, it is more instructive to consider the error plots. For a mesh of size 61 × 61 and SBP 42 in space, together with a sufficiently accurate SBP 84 in time, the errors at different times are presented in Figures 29–36.

To quantitatively compare the effects of the interfaces on accuracy, the error is integrated numerically over the spatial domain. The integrated errors for the two approaches when using SBP21 and SBP42 in space are shown in Figures 37 and 39. Additionally, we show the resulting CPU time for both methods in Figures 38 and 40. As can be seen, the gain from the new approach is twofold; we obtain lower error levels in the solution while spending less CPU time, even though the gain in CPU time for the 42 case is minimal. The reason for the large reduction in CPU time for the 21 case compared to the 42 case is presently not known to us. The convergence rates are of course identical.

(18)

The pressure distribution at t = 0.400 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 X -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Y 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 25. The standard multiblock.

The pressure distribution at t = 0.400

-0.2 0 0.2 0.4 0.6 0.8 1 1.2 X -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Y 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 26. The new approach.

-0.2 0 0.2 0.4 0.6 0.8 1 1.2 X -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Y 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 27. The standard multiblock.

-0.2 0 0.2 0.4 0.6 0.8 1 1.2 X -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Y 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 28. The new approach.

7.3. An application. As a final application, we consider a more complex geom-etry where no-penetration boundary conditions are imposed on the inner solid walls. At the outer boundaries, characteristic boundary conditions [3] with data from the manufactured solution

(53) U∞= [0, 0, 0, e−20((x−t)

2_+(y−t)2₎

]T are imposed.

A mesh of size 51 × 51 grid points in space and 201 nodes in time is constructed. Third and fifth order accurate SBP operators in space and time, respectively, are used. The pressure distribution and the velocity field at different times are shown in Figures 41–48. The resulting flow is tangential to the solid inner boundaries and the pressure pulse moves out of the domain as time passes.

8. Summary and conclusions. We have constructed a new combination of summation-by-parts operators which can be applied to a number of partial differential equations posed on non-simply connected spatial domains. To develop the theory, we considered a two-dimensional constant coefficient advection equation posed on a non-simply connected spatial domain and constructed an accurate and efficient scheme.

(19)

The error at t = 0.020 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y 1 2 3 4 5 6 7 8 9 10 #10-5

Fig. 29. The error, standard multiblock.

The error at t = 0.020 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y 1 2 3 4 5 6 7 8 9 10 #10-5

Fig. 30. The error, new approach. The error at t = 0.070 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y 1 2 3 4 5 6 7 8 9 10 #10-5

The error at t = 0.070 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y 1 2 3 4 5 6 7 8 9 10 #10-5

Fig. 32. The error, new approach. The error at t = 0.190 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y 1 2 3 4 5 6 7 8 9 10 #10-5

The error at t = 0.190 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y 1 2 3 4 5 6 7 8 9 10 #10-5

Fig. 34. The error, new approach.

Furthermore, we extended the new approach to a more complex non-simply connected geometries including curved domains. Although no proof of energy stability was obtained, correctly located spectra in combination with SBP in time indicates that the scheme is stable.

In the numerical experiments, we applied the new formulations to the linearized Euler equations. We showed that the new formulation is design order accurate by using the method of manufactured solutions. Additionally, we compared the error

(20)

The error at t = 0.213 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 x -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y 1 2 3 4 5 6 7 8 9 10 #10-5

Fig. 36. The error, new approach.

101 ₁₀2

Number of grid points in each space dimension, logarithmic 10-5

10-4

10-3

10-2 10-1

Energy norm of the errors, logarithmic

SBP21 in space, SBP84 in time

The new approach The standard technique

Fig. 37. Comparing the error levels.

101 ₁₀2

Number of grid points in each space dimension, logarithmic 10-1 100 101 102 103 104 105

CPU time [s], logarithmic

Fig. 38. Comparing the CPU times.

30 40 50 60 70 80 90

Number of grid points in each space dimension, logarithmic 10-6

10-5

10-4

10-3

Energy norm of the error, logarithmic

Fig. 39. Comparing the error levels.

30 40 50 60 70 80 90

Number of grid points in each space dimension, logarithmic 101

102

103

104

CPU time [s], logarithmic

Fig. 40. Comparing the CPU times.

levels and CPU time of the new approach with the standard multiblock technique. We conclude that the new method is more accurate and efficient compared with the standard multiblock technique. An application to a more complex geometry was also presented.

(21)