Simulation of wave propagation in discontinuous media and complex geometries

IT 10 056

Examensarbete 30 hp

Oktober 2010

Simulation of wave propagation

in discontinuous media and complex

geometries

Kristoffer Virta

Teknisk- naturvetenskaplig fakultet UTH-enheten Besöksadress: Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress: Box 536 751 21 Uppsala Telefon: 018 – 471 30 03 Telefax: 018 – 471 30 00 Hemsida: http://www.teknat.uu.se/student

Abstract

Simulation of wave propagation in discontinuous

media and complex geometries

Kristoffer Virta

A strictly stable high order finite difference (HOFD) scheme for the scalar wave equation with discontinuous coefficients corresponding to inhomogeneity in the underlying media is constructed using summation - by - parts (SBP) operators. The domain of the equation is allowed to be complex in the sense that it can not be smoothly mapped from the unit square. The inhomogeneous media and complex geometries are handled by treating the problem in a domain decomposition setting. For each sub-domain there is a smooth one to one mapping from the unit square and its underlying media is smooth. The sub-domains are then coupled together with a penalty method, the simultaneous approximation term (SAT) method. Finally the scheme is extended to be able to handle local grid refinement, a feature quite new in the context of the SBP - SAT methodology. Accuracy and convergence properties are verified in numerical computations in 2 dimensions.

Contents

1 Introduction 2

2 Definitions 3

2.1 The model problem . . . 4

2.2 Grid functions and related definitions . . . 4

2.3 Scalar products and norms . . . 6

2.4 Energy estimates and stability . . . 6

2.5 The summation - by - parts simultaneous approximation term method . . . 7

3 Analysis 9 3.1 Complex geometry . . . 9

3.2 The continuous case . . . 11

3.3 The semi - discrete case . . . 13

3.4 Non - conforming multi-block grids . . . 18

4 Numerical experiments 24 4.1 Case I - Cartesian coordinates . . . 24

4.2 Case II - Curvilinear coordinates and local grid refinement . . 26

5 Conclusions and future work 27

A The case of Cartesian sub - domains and piecewise constant

1

Introduction

Wave propagation problems arise in many applications, such as seismology, general relativity, acoustics and electromagnetics. In the general wave prop-agation problem the computational domain is large compared to the wave-length, so that the waves have to travel for long times (or correspondingly large distances). To efficiently resolve the solution during long times it is known that high order (at least third order) spatially accurate finite dif-ference schemes as well as a high order time marching method are ideally suited [1]. As the waves need to propagate during long times it is crucial to use a method that does not allow for any unphysical growth i.e. a strictly stable method. The media in which the waves travel is often discontinuous in practical applications. To model discontinuities without loosing accuracy and stability one must carefully consider how to treat the interfaces in the computational domain at which a discontinuity in the media occurs. Also, in general the region in which a solution is sought is complex in the sense that it can not be smoothly mapped from the unit square, which imposes difficulties if a finite difference scheme is to be used. Further, the underlying equations which are second order hyperbolic equations have historically with few exceptions been rewritten and solved as a system of first order equations. The obvious drawback with this approach is the doubling of the number of unknowns. It would therefor be beneficial to solve the equations on second order form.

The scalar wave equation written on second order form in discontinuous media have been approached successfully with a high order finite difference (HOFD) scheme on a Cartesian grid [9], in [7] a second order finite volume method is used to solve the acoustic wave equation in discontinuous media and arbitrary geometry. A second order FD scheme was used in [10] to solve the elastic wave equation in smooth media and non-complex geometries. This thesis aims at describing a strictly stable HOFD scheme to solve the scalar wave equation on second order form in discontinuous media allowing for complex geometries.

a smooth one to one mapping from the unit square and on which the media is smooth. As we shall see some thought needs to be invested in how to perform such a domain decomposition. Luckily the SBP - SAT methodology is well suited for such a task. Even though we can present our main result, a scheme that does meet the mentioned demands we shall observe that there is more to it. As waves progress over an interface of discontinuity the wave speed could be suddenly lowered yielding a more high frequent solution. To resolve the solution in the slower media while preserving a homogeneous accuracy we need to refine the grid locally. This is an obstacle which the SBP - SAT methodology has not been able to resolve, until recently. In [7] so called SBP preserving interpolation operators were presented, we shall show how to use these to resolve the issue of local grid refinement. To summarize, what we strive for is a method that:

• Is strictly stable.

• Solves the wave equation written on second order form. • Is explicit and of high order.

• Allows for discontinuous media. • Allows for complex geometries. • Allows for local grid refinement.

To limit the thesis in space and time we do not treat outer boundaries with any absorbing boundary conditions. However, the author acknowledge that this is of significant importance. We proceed as follows: in section 2 neces-sary definitions are provided, in section 3 we present our main result a scheme that resolves the first 5 items above together with a strict stability proof. As a corollary we can then construct with the newly developed SBP preserving interpolation operators a scheme that resolves the sixth item above. In sec-tion 4 we perform numerical experiments and in secsec-tion 5 we conclude and mention ideas for future work. Lastly an appendix is provided that hints on what future improvements may lead to.

2

Definitions

2.1

The model problem

In the following we will consider the 2 dimensional scalar wave equation. To allow for a piecewise continuous wave speed and geometries that can not be mapped smoothly from the unit square, from now on called complex geometries, we consider the equation in a domain decomposition setting. That is, the domain is a union of non-complex sub - domains on which the wave speed is smooth. Boundary conditions that are compatible at sub -domain corners are imposed on outer boundaries and as we shall see interface conditions at sub - domain interfaces are also needed for a stable problem. In other terms, we will consider the problem

utt = (bux)x+ (buy)_y, (x, y)∈ Ω

Liu = f, (x, y)∈ Ωj∩ Ωk (1)

Lbu = g, (x, y)∈ δΩ

Where Ω = ∪iΩi, b = bi > 0 (x, y),∈ Ωi, δΩ is the boundary of Ω, Li imposes sufficient interface conditions at the interfaces between sub - domains and

Lb imposes appropriate boundary conditions. An example of a two block scenario is illustrated in Figure 1.

Ω1

b = b1

Ω2

b = b2

γ

Figure 1: Ω = Ω1 ∪ Ω2. Special treatment at the interface Ω1 ∩ Ω2 = γ is

needed for a stable problem.

2.2

Grid functions and related definitions

We will work with rectangular sub - domains, either directly or by trans-forming the model problem on complex sub - domains to rectangular ones. Consider the rectangular domain Ω = {(x, y) | x1 ≤ x ≤ x2, y1 ≤ y ≤ y2}

Ω is discretized by introducing an equidistant N + 1 × M + 1 grid Ω =

6 y West -x South 6 u(0) North 6 u(k) 6 East u(N )

Figure 2: Rectangular domain

the following, to suppress notational clutter if not explicitly stated otherwise we let x2− x1 = y2 − y1 and N = M, so that hx = hy = h on all sub - do-mains. The value of a function u at a point (xi, yj) is denoted ui,j, this lets us introduce the vector u = [u(0)T_{, . . . , u}(N )T_]T _{where u}(k) _{= [u}

k,0, . . . , uk,N]T is the vector of values of u at xk along the y - direction, as illustrated in Figure 2. We will call u a grid function and throughout whenever a boldface notation is used grid functions are considered. With these conventions we can also treat grid functions as vectors and we define ui to be the element of

u at position i. Also, if we have defined a function f on a domain we mean

by a bold face f the corresponding grid function. If u(u) is a function(grid function) we let uE(uE), uW(uW), uN(uN) and uS(uS) denote the values of

u restricted to the east, west, north and south boundary (the vector of values

of u corresponding to the east, west, north and south boundary). Let v be a vector by d(v) we mean the diagonal matrix with the elements of v on the diagonal. We will make use of the Kroenecker product of two matrices A and

B of size p× q and r × s, C = A⊗ B =    a1,1B . . . a1,qB .. . . .. ... ap,1B . . . ap,qB    .

The following properties of the Kroenecker product will be used

(A⊗ B)(C ⊗ D) = (AC) ⊗ (BD), whenever AC and BD are defined (A⊗ B)−1 = (A−1⊗ B−1), whenever A−1 and B−1 exist

(A⊗ B)T = AT ⊗ BT.

vector of the same size as a grid function on the domain. The insertion should be at the corresponding boundary positions and the resulting vector should have zeros elsewhere. For this purpose the following matrices are useful

iE = eR⊗ I, iW = eL⊗ I

iN = I ⊗ eU, iS = I ⊗ eD (2) Where eR, eL, eU, eD are the following vectors of length N

eR = eU =    0 .. . 1    , eL = eD =    1 .. . 0    .

That is, iEuE gives a vector of the same length as u but with non - zero elements only at the positions corresponding to the east boundary.

2.3

Scalar products and norms

For functions u, v ∈ L2(Ω) we will use the scalar product (u, v)w =

∫

Ω

uvwdx, w(x) > 0

With the corresponding norm

∥u∥2

w = (u, u)w

when u, v are grid functions and H a symmetric positive definite matrix we will use the scalar product

(u, v)H = uTHv With the corresponding norm

∥u∥2

H = (u, u)H

2.4

Energy estimates and stability

Consider the one dimensional scalar wave equation

utt = (bux)x, b > 0, t≥ 0, x ∈ [x1, x2]. (3)

The kinetic energy of (3) is u2_t/2 and the potential energy is bu2_x/2 why the

Multiplying (3) by ut and integrating by parts (the energy method) leads to

d

dtE(t) = 2utbux|

x2

x1 (5)

Hence, the behaviour of the total energy over time is governed by the bound-ary conditions at x1 and x2. Taking for example homogeneous Neumann

conditions yields

d

dtE(t) = 0, ⇒ E(t) ≤ C, C constant. (6)

We call the inequality (5) an energy estimate. We make the analogous defini-tions for the model problem. That is, we define an energy of the equation, a non - negative integral containing a interpretation of the kinetic energy and the potential energy and say that we have an energy estimate if the time be-haviour of the energy is governed by the outer boundary conditions. For any numerical scheme that solves (1) to be useful the scheme must at least nearly satisfy an analogous discrete energy estimate. We shall construct a semi -discretization of (1) and define an discrete energy analogous to the continuous one and prove that it mimics the continuous energy estimate exactly. Hence, rendering any unphysical growth in the numerical solution impossible. Such a scheme is called strictly stable and the method of mimicking the continuous energy estimate is called the discrete energy method. In general there are many different definitions of stability of initial boundary value problems see for example chapters 9 - 12 in [2].

2.5

The summation - by - parts simultaneous

approxi-mation term method

Summation by parts (SBP) operators are discrete operators approximating spatial derivatives that mimic the continuous integration by parts property. In particular, they allow for construction of semi - discrete schemes that satisfies a discrete energy estimate mimicking (5). High order SBP operators for the first derivative has previously been developed in for example [3]. An SBP operator D1 approximating the first derivative has the following

properties

D1 = H−1Q, H = HT > 0

SAT method imposes boundary conditions weakly by adding a forcing term to the spatial discretization which is proportional to the difference between the discrete and imposed value at the boundary. Consider for example the problem (3) with Neumann boundary conditions ux(t, 0) = f1, ux(t, 1) = f2.

If u is the one - dimensional grid function corresponding to u of (3), then the semi - discretization utilising the SBP - SAT technique can be written

utt = D1BD1u + τ1H−1eL((D1u)1− f11) + τ2H−1eR((D1u)N − f2N), (8)

where B = d(b). To illustrate the power of the SBP - SAT methodology we show that we can obtain a discrete energy estimate that mimics (5). Multiplying (8) by uT_tH, adding the transpose (the discrete energy method)

and using the properties of (7) yields

uT_tHutt+ uTttHut=−2(D1ut)THB(D1u)

+ 2 (utN(BD1u)N + ut0(BD1u)0)

+ 2τ1ut0(BD1u)0− 2τ1ut0f10

+ 2τ2utN(BD1u)N − 2τ2utNf2N

Now, taking τ1 =−1, τ2 =−1 gives

d dt ( ∥ut∥H + (D1u)THBD1u ) = 2(utNf2N − ut0f10)

Hence, if we define the discrete energy as

E(t) =∥ut∥2H +∥D1u∥2HB

we get an discrete energy estimate mimicking (5). Taking homogeneous boundary conditions, f1 = f2 = 0 gives us the discrete analogue of (6).

As is seen in an example in [8] it is possible to initially choose quite arbitrary boundary values and still get a solution with boundary values that converges to the true ones. This property of the SAT method to make the boundary values converge to the imposed boundary value regardless of the initial value has given it the folklore nickname "rubber band technique". The parameters arising in the SAT method are called penalty parameters and it is a more or less arduous task to determine their appropriate values. The matrix H can be both a full or a diagonal matrix. In the following we will use only operators

D1 with the corresponding H diagonal. See [3] for details of full norm SBP

operators. Let (x1, x2) be a Cartesian coordinate system, to approximate the

partial derivatives _∂x∂

1,

∂

∂x2, we will use the operators

we will also use the following matrices

Hx1 = H ⊗ I, Hx2 = I⊗ H.

We note how

Hx1x2 = Hx1Hx2

defines a norm for two dimensional grid functions.

3

Analysis

This is the main part of the work. We begin by describing how to transform the model problem to a problem that can be solved on Cartesian sub - do-mains. Then we determine sufficient conditions that needs to be imposed at sub - domain interfaces to obtain a stable problem. After that we utilise the SBP - SAT methodology to describe a semi - discretization of the continuous case and prove that it is strictly stable by mimicking the continuous energy estimate. Finally, as a corollary we can present a scheme that enables local grid refinement.

3.1

Complex geometry

To allow for complex geometries we will proceed with curvilinear coordinates. We begin by transforming the equation of (1) on piecewise curvilinear do-mains to an equation on piecewise rectangular sub - dodo-mains. Then we derive sufficient interface conditions for an suitable energy estimate in the continu-ous case and then mimic these in a semi - discretization, for which we show stability by deriving a discrete energy estimate that mimics the continuous one. We develop the theory for a domain consisting of two sub - domains coupled at the east respective west boundaries, since as we shall see interface conditions are local. Meaning, that if we prove strict stability for the case of two sub - domains, stability follows for an arbitrary number of sub - domains. Transforming the equation

Suppose that there is a a one-to-one mapping

(x, y) = (x(ξ, η), y(ξ, η)), (ξ, η)∈ [ξ1, ξ1+ 1]× [η1, η1+ 1]

from a rectangle of side one, the computational domain, to a sub - domain where we wish to solve the model problem, the physical domain, figure 3. The Jacobian (J ) of the mapping is

-y 6 x




 ξ  * η  x = x(ξ, η) y = y(ξ, η) -η 6 ξ

Figure 3: The mapping between the Cartesian (right) and curvilinear (left) coordinates

The Jacobian is everywhere non - zero since the mapping is one-to-one. By the chain rule we have

{ uξ = uxxξ+ uyyξ uη = uxxη+ uyyη ⇔ { ux = uξyη−uηyξ J uy = −uξ xη+uηxξ J ⇔ { ux = (uyη)ξ−(uyξ)η J uy = (−uxη)ξ+(uxξ)η J Using these relations to calculate the second derivatives we get

   uxx = _J1 ( (uξyη−uηyξ) J yη ) ξ− 1 J ( (uξyη−uηyξ) J yξ ) η uyy = −1_J ( (uηxξ−uξxη) J xη ) ξ+ 1 J ( (uηxξ−uξxη) J xξ ) η Summing (bux)x+ (buy)y and rearranging terms (1) transforms into

J utt = (bF u)ξ+ (bGu)η. (9) Where F = α∂ξ+ β∂η, G = γ∂η + β∂ξ α = x 2 η + yη2 J , β =− xξxη + yξyη J , γ = x2_ξ+ y_ξ2 J .

That is, if we solve the equation (9) on the rectangle [ξ1, ξ1+ 1]× [η1, η1+ 1]

y x -6 utt= b(uxx+ uyy) Ω1 b = b1 x1= x1(ξ, η) y1= y1(ξ, η) Ω2 b = b2 x2= x2(ξ, η) y2= y2(ξ, η) J utt= (bF u)x+ (bGu)x Ω′₁ F = F1, G = G1 u = v Ω′₂ F = F2, G = G2 u = w η ξ -6

Figure 4: The left physical domain and the right computational domain Ωi, there is a one - to - one mapping from a rectangle, Ω′i such that adjacent rectangles are mapped to adjacent sub - domains. Thus, on rectangle i we solve Jiutt = (biFiu)ξ+(biGiu)η, where Ji, Fi, Gi are defined by the mapping from Ω′_i to Ωi.

3.2

The continuous case

Let Ω1 be the left sub - domain and Ω2 the right sub - domain. We assume

that there is a one - to - one mapping from the rectangle Ω′₁ = [−1, 0]×[0, 1] to Ω1 that defines α1, β1, γ1, and thus F1 and G1. And a one - to - one

mapping form the rectangle Ω′₂ = [0, 1]× [0, 1] to Ω2 that defines α2, β2,

γ2, and thus F2 and G2. We shall solve the transformed equation (9) on the

computational domain Ω′ = Ω′₁ ∪ Ω′₂ to resolve the model problem on the physical domain Ω = Ω1 ∪ Ω2. We will denote by v the solution on Ω′1 and

method applied to the equation (9) gives d dt (∥vt∥J1 +∥wt∥J2) = ∫ Ω′ ut((bF u)ξ+ (bGu)η) dξdη = ∫ 1 0 vtb1F1v|0₋₁dη − ∫ Ω′₁ vtξb1F1vdξdη + ∫ 0 −1 vtb1G1v|10dξ− ∫ Ω′₁ vtηb1G1vdξdη + ∫ 1 0 wtb2F2w|10dη − ∫ Ω′₂ wtξb2F2wdξdη + ∫ 1 0 vtb2G2w|10dξ− ∫ Ω′₂ wtηb2G2wdξdη =− ∫ Ω′1 vtξb1(α1vξ+ β1vη) + vtηb1(γ1vη+ β1vξ) dξdη − ∫ Ω′₂ wtξb2(α2wξ+ β2wη) + wtηb2(γ2wη + β2wξ) dξdη + (wtE, (F2w)E)b_2E − (vtW, (F1v)W)b_1W + (utN, (Gu)N)bN − (utS, (Gu)S)bS + ∫ 1 0 vtE(b1F1v)E− wtW(b2F2w)Wdη. Now, vtξb1(α1vξ+ β1vη) + vtηb1(γ1vη+ β1vξ) = b1 d dt ( α1vξ2+ 2β1vξvη + γ1vη2 ) = b1 d dt ([ vξ vη ] [α1 β1 β1 γ1 ] [ vξ vη ])

Further, since α1 > 0 and α1γ1−β12 = (x1ξy1η−x1ηy1ξ)2 = J12 > 0, the matrix

[

α1 β1

β1 γ1

]

is positive definite. The analogous result holds for the integral containing only w. Thus, we can define the energy as

We have d dtE(t) = (wtE, (F2w)E)b2E − (vtW, (F1v)W)b1W + (utN, (Gu)N)bN − (utS, (Gu)S)bS + ∫ 1 0 vtE(b1F1v)E − wtW(b2F2w)Wdη

We see that to obtain an energy estimate that is governed by the bound-ary conditions we need the following conditions vtE = wtW, (b1F1v)E = (b2F2w)W. However, vtE = wtW as well as the equality of time derivatives of any order at the interface is necessary for vE = wE. That is, if we impose the interface conditions

vE = wW, (b1F1v)E = (b2F2w)W (10) then we get the energy estimate

d

dtE(t) = (wtE, (F2w)E)b2E − (vtW, (F1v)W)b1W (11)

+(utN, (Gu)N)bN − (utS, (Gu)S)bS

3.3

The semi - discrete case

We will now describe the semi - discretization of (9) with interface condi-tions that are consistent with (10). We will impose homogeneous Neumann boundary conditions

F u = 0, ξ, η∈ Ω′_E ∪ Ω′_W, Gu = 0, ξ, η∈ Ω′_N ∪ Ω′_S

on the outer boundaries. This yields the continuous energy estimate

d

dtE(t) = 0 (12)

To prove strict stability a discrete energy estimate mimicking (12) for the semi - discretisation will be derived. The discretized F1,2 and G1,2 look like

F1,2 = ( d(α1,2)D1ξ + d(β1,2)D1η ) G1,2 = ( d(γ1,2)D1η+ d(β1,2)D1ξ ) .

The semi - discretization of the interface conditions (10) can be written

where B1 = d(b1) and B2 = d(b2). In the continuous case the continuity of

time derivatives of any order across sub-domain interfaces follows from the continuity of the solution across interfaces. In the semi-discretization we can impose the continuity of time derivatives explicitly in order to investigate if any gains are to be made in terms of accuracy and convergence by doing so. The semi-discretizatin of the explicit imposition of the continuity of first and second time derivatives can be written

I3 = vtE− wtW = 0, I4 = vttE − wttW = 0.

To impose the continuity of the solution over the interface with the SAT method we introduce the penalty matrix τ through

τ =    τ1 . . . 0 .. . . .. ... 0 . . . τN    The following matrices are also needed

R =    0 . . . 0 .. . . .. ... 0 . . . 1    , L =    1 . . . 0 .. . . .. ... 0 . . . 0   

The semi - discrete scheme corresponding to the continuous case using SBP operators to approximate the spatial derivatives and the SAT method to impose outer boundary and interface conditions can now be written:

J1vtt = DξB1F1v + DηB1G1v J2wtt = DξB2F2w + DηB2G2w + H_ξ−1(R⊗ τ)iEI1 − Hξ−1(L⊗ τ)iWI1 +1 2H −1 ξη (B1F1)THηiEI1 − 1 2H −1 ξη (B2F2)THξiWI1 +1 2H −1 ξ iEI2 − 1 2H −1 ξ iWI2 (13) + σH_ξ−1iEI3 − σH_ξ−1iWI3 + ρH_x−1iEI4 − ρHx−1iWI4 + H_ξ−1iWB1F1v − Hξ−1iEB2F2w − H−1 η iNB1G1v − Hη−1iNB2G2w + H_η−1iSB1G1v + Hη−1iSB2G2w

the continuity of first and second time derivatives across the sub - domain interfaces. We will below refer to these terms as the extra conditions.

We are now ready to prove the main result of the present study: Theorem 1 Let h1 = H1,1, hN = HN,N if τi ≤ −

b_1Eiα_1Ei

4hN −

b_2Wiα_2Wi

4h1 , σ ≤ 0

and ρ≤ 0, then the scheme (13) is strictly stable.

Proof:

The discrete energy method applied to the scheme (13) gives

d dt ( ∥vt∥J1Hξη +∥wt∥J2Hξη + Y T 1 A1Y1+ Y2TA2Y2 ) = 2σZTCZ + d dt ( XTM X + ρZTCZ) Where Y1 = [ D1ξv D1ηv ] , Y2 = [ D1ξw D1ηw ] A1 = ¯HξηB¯1C1, A2 = ¯HξηB¯2C2, ¯ B1 = [ B1 0 0 B1 ] , ¯B2 = [ B2 0 0 B2 ] , C1 = [ d(α1) d(β1) d(β1) d(γ1) ] , C2 = [ d(α2) d(β2) d(β2) d(γ2) ] , ¯Hξη = [ Hξη 0 0 Hξη ] , X =     vE wW (F1v)E (F2w)W     , M =     τ H −τH 1/2H 1/2H −τH τ H −1/2H −1/2H 1/2H −1/2H 0 0 1/2H −1/2H 0 0     , Z = [ vtE wtE ] and C = [ 1 −1 −1 1 ] ⊗ H.

C is positive semi-definite why taking σ ≤ 0 and ρ ≤ 0 guarantees that the

terms containing Z are always non-positive. We prove the positive definite-ness of A1 and A2, in order to be able to define a sensible discrete energy.

Let z be any vector of length 2N2_{, then}

To prove a strict inequality suppose there is a z non - zero such that zTC1z =

0, then for all i

(x_{η i}zi− x_{ξ i}zN2_+i)2 = 0 and (y_{η i}z_i− y_{ξ i}z_N2_+i)2 = 0.

That is for all i,

x_{η i} x_{ξ i} = zN2_+i zi = yη i y_{ξ i} ⇔ 0 = x_{ξ i}y_{η i}− x_{η i}y_{ξ i} = J1i

which is a contradiction. Further, ¯H ¯B1C1 =

√ ¯ H ¯B1C1 √ ¯ H ¯B1 the positive

definiteness of A1 follows and similarly for A2. Now, if M is negative

semi-definite we could define the discrete energy as

E(t) =∥vt∥J1Hξη +∥wt∥J2Hξη + Y

T

1 A1Y1+ Y2TA2Y2

−ρZT_{CZ − X}Y_{M X.}

However, M is indefinite. Let c1 and c2 be diagonal matrices of size N× N

and ˜ M =     τ H −τH 1/2H 1/2H −τH τ H −1/2H −1/2H 1/2H −1/2H −c1H 0 1/2H −1/2H 0 −c2H     . We have XTM X = X˜ TM X − Y₁TA₁′Y1− Y2TA′2Y2 Where A′₁ = ¯HηB¯1 2 C₁′, A′₂ = ¯HηB¯2 2 C₂′ C₁′ = [ (R⊗ c1)α12 (R⊗ c1)α1β1 (R⊗ c1)β1α1 (R⊗ c1)β12 ] , C₂′ = [ (L⊗ c2)α22 (L⊗ c2)α2β2 (L⊗ c2)β2α2 (L⊗ c2)β22 ] ¯ Hξ = [ Hξ 0 0 Hξ ] , ¯Hη = [ Hη 0 0 Hη ] Introduce ˜ A1 = A1− A′1, ˜A2 = A2− A′2

If we can choose c1 and c2 such that ˜A1 and ˜A2 are at least positive semi

-definite, then if we can choose τ to make ˜M negative semi - definite we can

define the discrete energy to be ˜

E(t) =∥vt∥J1Hξη +∥wt∥J2Hξη + Y

T

1 A˜1Y1+ Y2TA˜2Y2

and get the energy estimate

d

dtE(t) = 2σZ˜

T_CZ

which mimics (12) with if σ < 0 some additional damping. Now ˜ A1 = ¯HηB¯1C˜1, ˜C1 = (_¯ HξC1− ¯B1C1′ ) .

If ˜C1 is positive semi - definite then so is ˜A1. We consider the Cholesky

factorisation of ˜C1

˜

C1 = GGT,

where G is lower triangular. Since ˜A1 is tridiagonal it is particularly simple

to write down the explicit Cholesky factorisation of ˜C1. The elements of G

affected by c1 are:        Gj,j = (hNα1Ei− c1ib1Eiα 2 1Ei)1/2 G_N2_+j,j =

hNβ1Ei−c1ib_1Eiα1Eiβ1Ei

(hNα1Ei−c1ib_1Eiα2_1Ei)1/2

G_N2_+j,N2_+j = ( hNγ1Ei− c1ib1Eiβ 2 1 Ei− ( _h

Nβ1Ei−c1ib_1Eiα1Eiβ1j

(hNα1Ei−c1ib_1Eiα2_1Ei)1/2

))1/2

,

where

j = N2− N + 1 . . . N2, i = j− N2+ N.

We know that G having real non - negative diagonal elements is equivalent to ˜C1 = GGT being positive semi - definite. See for example [6]. Hence, we

want to choose c1 such that the diagonal elements of G are real and non

-negative. Let now c1i→ _b hN

1Eiα1Ei, then

Gj,j → 0 GN2_+j,j = β_1Ei α1 1/2 j hN − c1ib1Eiα1Ei (hN − c1ib1Eiα1Ei) 1/2 → 0 GN2_+j,N2_+j → ( hN(γ1Ei− β12Ei α1Ei ) )1/2 .

GN2_+j,N2_+j is real and positive since α_1Eiγ_1i−β₁2_Ei = J_1i > 0. If c_1iis at most

equal to hN

β1Eiα1Ei, then ˜A1 is positive semi - definite. In an analogous way we

determine c2 such that ˜A2 is at least positive semi - definite. Investigation

with the Matlab symbolic toolbox shows that taking

renders ˜M negative semi - definite 

We introduce the penalty parameter ΓCu through

τ = ΓCu ( −d(b1E)d(α1E) 4hN −d(b2W)d(α2W) 4h1 ) Time integration

The semi - discretization (13) can be written on the form

utt = F u + Gut, if σ ̸= 0, ρ ̸= 0 (14) Where for the case σ = ρ = 0 G reduces to the zero matrix. To solve (14) we will use the classic Runge - Kutta 4th order explicit method (RK4). For this we put v = ut and solve with RK4 the problem

[ u v ] t = [ 0 I F G ] [ u v ] , [ u(x, y, 0) v(x, y, 0) ] = [ u0 v0 ] , (15) where I is the N2 × N2 identity matrix. The stability proven above can be verified numerically by plotting the eigenvalues of the system matrix in (15). For stability the eigenvalues must be in the stability region of RK4 i.e. have non - positive real part, the magnitude of the spectral radius is handled by adjusting the time step. In Figure 5 we have chosen N = 26,

b1 = 1, b2 = 1/2. To show the importance of choosing the right parameters

figure 5(b) and 5(d) displays how an slightly erroneous choice of parameters destroys the stability, note the difference in scale on the x - axis.

3.4

Non - conforming multi-block grids

−5 0 5 x 10−14 −40 −30 −20 −10 0 10 20 30 40 Real Imag (a) ΓCa= 1, σ = ρ = 0 −1.5 −1 −0.5 0 0.5 1 1.5 x 10−5 −40 −30 −20 −10 0 10 20 30 40 Real Imag

(b) σ = ρ = 0. Multiplying one diagonal element of c1by 0.99 renders the scheme

un-stable −0.4 −0.35 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 −40 −30 −20 −10 0 10 20 30 40 Real (c) ΓCu= 1, σ =−0.1, ρ = −0.1 −2 0 2 4 6 8 10 12 14 x 10−3 −40 −30 −20 −10 0 10 20 30 40 Real Imag (d) ΓCu = 1, choosing σ = 0.001, ρ = 0.001

destroys the stability

Figure 5: The left column shows the spectrum of the system matrix of (15) with a stable parameter choice, the right column with parameters that violate Theorem 1

domains arises also in other problems and has limited the use of the SPB -SAT methodology in for example adaptive solvers. Recently Mattsson et.al constructed in [7] interpolation operators that if used properly does not de-stroy stability, SBP preserving interpolation operators. We shall show how to use these for our model problem. We will consider a domain consisting of two sub - domains with a doubling of the resolution of the grid going from the left, coarse grid, to the right, fine grid, sub - domain, (see Figure 8). Let

the coarse respective fine grid, we introduce the following notation: IC, the NC× NC identity IF, the NF × NF identity D1ξC = D1C ⊗ IC, D1ηC = IC ⊗ D1C D1ξF = D1F ⊗ IF, D1ηF = IF ⊗ D1F iEC = eC⊗ IC, iW F = eF ⊗ IF

Where D1C = HC−1QC and D1F = HF−1QF are 1 - D SBP operators approx-imating the first derivative on grid with NC respective NF grid points, eE and eF are the following vectors of length NC respective NF,

eC =    0 .. . 1    , eF =    1 .. . 0    .

Now, we shall use the interpolation operators IC2F, that elongate a grid function from the coarse grid to the fine grid, and IF 2C, that restricts a grid function from a fine grid to a coarse grid. The interpolation operators have (see [7] for details) the properties

HFIC2F = IF 2CT HC (16)

HC(IC − HCIF 2CIC2F)≥ 0, HF(IF − HFIC2FIF 2C)≥ 0.

The properties (16) is as we shall see crucial for stability. We denote by v the solution on the left coarse grid and by w the solution on the right fine grid. Further, we assume one - to one mappings from the rectangles [−1, 0]×[0, 1] and [0, 1]× [0, 1] to the left and right grid that define JC, JF, FC, FF, GC

and GF. Using the interpolation operators the semi - discretization of the

interface conditions (10) becomes

I1C = vE − IF 2CwW = 0, I1F = wW − IC2FvE = 0, and

I2C = (BCFCv)E − IF 2C(BFFFw)W = 0,

I2F = (BFFFw)W − IC2F(BCFCv)E = 0. The extra conditions can be written

I3C = vtE− IF 2CwtW = 0, I3F = wtW − IC2FvtE = 0,

We can now write the semi - discretization as: JCvtt = D1ξCBCFCv + D1ηCBCGCv JFwtt = D1ξFBFFFw + D1ηFBFGFw + τ H_Cξ−1iECI1C − τH_{F ξ}−1iW FI1F +1 2H −1 Cξη(BCFC) T_H CηiECI1C − 1 2H −1 F ξη(BFFF) T_H F ξiW FI1F +1 2H −1 CξiECI2C − 1 2H −1 F ξiW FI2F + σH_Cξ−1iECI3C − σHF ξ−1iW FI3F + ρH_Cξ−1iECI4C − ρHF ξ−1iW FI4F (17) + H_Cξ−1iW CBCFCv − HF ξ−1iEFBFFFw − H−1 CηiN CBCGCv − HF η−1iN FBFGFw + H_Cη−1iSCBCGCv + H_{F η}−1iSFBFGFw

The proof of stability is similar to the proof of theorem 1, we have: Corollary 1 Let hF = HF1,1, hC = HC_{NC ,NC} if τ ≤ − 1 4cC − 1 4cF where cC = min 1≤i≤NC { hN bCiαCi}, CF = _1≤i≤Nmin_F{ h1

bFiαFi}, σ ≤ 0 and ρ ≤ 0 then the scheme

(17) is strictly stable.

Proof:

To improve readability we consider the case σ = 0 and ρ = 0 the stability of case σ < 0 and ρ < 0 is proven analogously to the proof in Theorem 1 using the properties (16) in a way similar to what is done below. The discrete energy method yields

d dt ( ∥vt∥JC +∥wt∥JF + Y T CACYC + YFTAFYF ) = d dt ( τ v_ETHCvE − 2vTEHCIF 2CwW + τ wWT HFwW ) + d dt ( (FCv)TEHCvE − (FFw)TWHFwW ) − (FCvt)TEHCIF 2CwW − wTtWHFIC2F(FCv)E + (FFwt)TWHFIC2FvE + vtEHCIF 2C(FFw)W We now use the first of the properties (16) to get

wT_tWHFIC2F(FCv)E = (IF 2CwtW)THC(FCv)E and

Thus, we have d dt ( ∥vt∥JC +∥wt∥JF + Y T CACYC + YFTAFYF − XTM X ) = 0 Where YC = [ D1Cξv D1Cηv ] , YF = [ D1F ξw D1F ηw ] , AC = ¯HCξηB¯CCC, AC = ¯HF ξηB¯FCF, ¯ BC = [ BC 0 0 BC ] , ¯BF = [ BF 0 0 BF ] , CC = [ d(αC) d(βC) d(βC) d(γC) ] , CF = [ d(αF) d(βF) d(βF) d(γF) ] , ¯ HCξη = [ HCξη 0 0 HCξη ] , ¯HF ξη = [ HF ξη 0 0 HF ξη ] , X =     vE IF 2CwW (BCFCv)E IF 2C(BCFFw)W     , M =     τ −τ 1/2 1/2 −τ τ −1/2 −1/2 1/2 −1/2 0 0 1/2 −1/2 0 0     ⊗ HC.

We are in a situation very similar to before. We define

˜ M =     τ −τ 1/2 1/2 −τ τ −1/2 −1/2 1/2 −1/2 −cC 0 1/2 −1/2 0 −cF     ⊗ HC We have XTM X = X˜ TM X − cCYCA′CYC − cF(IF 2C(BFFFw)W)THCIF 2C(BFFFw)Ww)W where A′_C = ¯HηB¯C2CC′ , C_C′ = [ (RC⊗ IC)α2C (RC⊗ IC)αCβC (RC ⊗ IC)βαC (RC⊗ IC)β2C ] , ¯Hη = [ HCη 0 0 HCη ] .

min

1≤i≤NC

{ hN

bCiαCi} must be taken for cC. Further, we have that

− cF(IF 2C(BFFFw)W)THCIF 2C(BFFFw)W =−cF(BFFFw)TWHFIC2FIF 2C(BFFFw)W = cF(BFFFw)TWHF(IF − IC2FIF 2C)(BFFFw)W − cF(BFFFw)TWHF(BFFFw)W. Now −cF(BFFFw)TWHF(BFFFw)W =−cFYFA′FYF Where A′_F = ¯HηB¯F2CF′ , C_F′ = [ (LF ⊗ IF)α2F (LF ⊗ IF)αFβF (LF ⊗ IF)βαF (LF ⊗ IF)βF2 ]

and LF is the NF×NF matrix with LF1,1 = 1, LFi,j = 0, (i, j)̸= (1, 1). Again

as in theorem 1 we define ˜AF similarly to ˜A2 and determine cF similarly to c2

with the difference that now cF is a scalar we must take cF = min

1≤i≤NF

{ h1 bFiαFi}

why we can define the discrete energy as ˜ E(t) =∥vt∥JC +∥wt∥JF + bCY T CA˜CYC+ bFYFTA˜FYF −XT_{M X + c˜} F(BFFFw)TWHF(IF − IC2FIF 2C)(BFFF)W.

Here we use the second property of (16) and the fact that cF > 0 to see that the energy is always non - negative. We get the energy estimate

d dt

˜

E(t) = 0 △

We introduce the penalty parameter ΓCF through

τ = ΓCF(− 1 4cC − 1 4cF )

Remark: Why not introduce a vector of penalty parameters as in theorem 2, it gives a sharper boundary for the penalty parameters resulting in lesser stiffness in the resulting system of ordinary differential equations and a more elegant theory. If we look at the following term arising in the proof for stability:

−2wT

tWτFHFIC2FvE

4

Numerical experiments

In this section we will perform numerical experiments to verify the theory developed above. In a first experiment we investigate the influence of the extra conditions on accuracy and convergence, for this we consider the model problem on a rectangular domain. In a second experiment we consider the model problem on a complex domain where we also apply the technique of local grid refinement. In both experiments comparisons are made between the numerical solution uapprox and the analytic solution projected onto the

grid, uexact. We measure the discrete L2 error

e(N ) =∥uapprox− uexact∥h.

where uapprox is approximated on a grid of size corresponding to N and uexact

is resolved on the same grid. We measure the rate of convergence q as

q = log10 e(N1) e(N2) log₁₀ ( N2 1 N2 2 )1/d.

4.1

Case I - Cartesian coordinates

To investigate if any gains are to be made by the explicit imposition of the continuity of the first and second time derivative across a discontinuity, the extra conditions, we solve the model problem on rectangular sub - domains with a discontinuous coefficient. We chose the domain [−1, 0] × [0, 1] ∪ [0, 1]× [0, 1] with a block interface at x = 0. Further b1 = 1/2, b2 = 9/34

and homogeneous Neumann boundary conditions are imposed. The model problem then has the solution

u(x, y, t) =

{

cos(3πt) cos(3πx) cos(3πy), (x, y)∈ [−1, 0] × [0, 1] cos(3πt) cos(5πx) cos(3πy), (x, y)∈ [0, 1] × [0, 1]

We use RK4 for time integration with a time step of 10−4 and compute until

t = 0.9. SBP operators with a fourth and sixth order interior stencil are

SBP operator for the first derivative 6th order accurate in the interior and 4th order accurate at the boundaries two times to approximate the second derivative an 4th order accurate approximation for the second derivative should result.

From the results of the experiment we note how one should chose ΓCu to be slightly above the minimal allowed value 1 to get optimal accuracy and order of convergence. Further, we see that using the 4th order stencil gives a rate of convergence higher than the theoretical expected value. Finally, we also see that, unfortunately, no gains are to be made from the explicit imposition of the first and second time derivative across the interface.

4th order interior stencil

ΓCu = 1 ΓCu = 1.2 ΓCu = 5

N log₁₀e(N ) q log₁₀e(N ) q log₁₀e(N ) q

26 -1.52 - -1.52 - -1.54

-51 -2.69 4.00 -2.81 4.39 -2.83 4.42

101 -3.69 3.35 -4.01 4.04 -4.05 4.10

201 -4.61 3.08 -5.20 4.01 -5.26 4.03

6th order interior stencil

ΓCu = 1 ΓCu = 1.2 ΓCu = 5

N log₁₀e(N ) q log₁₀e(N ) q log₁₀e(N ) q

26 -1.54 - -1.65 - -1.63

-51 -2.86 4.50 -2.88 4.19 -2.88 4.27

101 -4.17 4.44 -4.12 4.19 -4.11 4.14

201 -5.51 4.46 -5.43 4.37 -5.43 4.41

Table 1: Rate of convergence and accuracy for different values of ΓCawithout the extra conditions

ΓCu = 2, σ =−4, ρ = −4

4th order interior stencil 6th order interior stencil N log₁₀e(N ) q log₁₀e(N ) q

26 -1.54 - -1.63

-51 -2.83 4.42 -2.88 4.28

101 -4.05 4.11 -4.11 4.14

201 -5.26 4.04 -5.43 4.41

4.2

Case II - Curvilinear coordinates and local grid

re-finement

To test the accuracy and rate of convergence on a problem in curvilinear coordinates using local grid refinement we choose an analytic solution which is defined on the whole ofR2_{, with discontinuities in the coefficient along the}

lines x− y = n, n ∈ Z.

u(x, y, t) =

{

cos(3πt) cos(2π(x− y)), (x, y) ∈ S2n+1

cos(3πt) cos(4π(x− y)), (x, y) ∈ S2n

Where S2n+1is the strip{(x, y)|−∞ < x < ∞, x−(2n+1) ≤ y ≤ x−2(n+1)}

and S2n is the strip {(x, y)| − ∞ < x < ∞, x − (2n + 1) ≤ y ≤ x −

2n}. To solve this problem numerically we proceed with periodic boundary conditions. We use the domain depicted in Figure 6. On S2n the solution is

more high frequent why we apply a local grid refinement on the right sub -domain. As we saw in the above experiment imposing the extra conditions did not yield any improvement, why we perform tests with σ = ρ = 0 in this experiment. We use RK4 with a time step of 10−4 for time integration and we compute until t = 0.8. The order of accuracy and rate of convergence of using schemes of 4th and 6th order in the interior is displayed in Table 3. One can show that approximating _dxd(f_dxd) when f is variable, which we do in this experiment due to the transformation to curvilinear coordinates, with an SBP operator approximating the first derivative 4th order accurate in the interior and 2nd order accurate at the boundaries results in an 4th order accurate approximation of _dxd(f_dxd). Similarly, using an SBP operator 6th order accurate in the interior and 4th order accurate at the boundaries results in an 5th order accurate approximation of _dxd(f_dxd). The entries for the rate of convergence is deviant in the row corresponding to NC = 51 probably a result of the numerical solution not being sufficiently resolved for NC = 51.

ΓCF = 2

4th order interior stencil 6th order interior stencil

NC log10er(NC) q(NC) log10er(NC) q(NC)

26 -1.54 - -1.68

-51 -2.65 3.79 -2.55 2.96

101 -3.84 4.03 -4.12 5.29

151 -4.54 3.96 -4.99 5.02

−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6

Figure 6: We impose periodic boundary conditions on the outer boundary.

5

Conclusions and future work

In this thesis we have described a HOFD scheme for solving the scalar wave equation in discontinuous media and complex geometries. We also proved the strict stability of the scheme. We noted that in order to maintain uni-form accuracy we needed the possibility to peruni-form local grid refinement. We showed how this could be done with the use of SBP preserving inter-polation operators without loosing strict stability or accuracy. Numerical experiments verified the theory. In an appendix we considered the special case of piecewise constant media and a rectangular geometry and saw that for this case we could achieve considerably better results with the use of narrow diagonal SBP operators. From this result a conjecture is made that if the corresponding narrow diagonal SBP operators can be constructed for general variable coefficients, which is the case in curvilinear coordinates and a general media, then the corresponding improvement is to be made. The investigation of this conjecture is a natural continuation of the present work and will be handled in an upcoming study. Other continuations could in-volve extending the scheme to resolve 3 dimensional geometries or solve for the elastic, electromagnetic, etc wave equations.

Acknowledgements

Appendix

A

The case of Cartesian sub - domains and

piecewise constant coefficients

We shall now consider the special case in which all sub - domains are Carte-sian blocks and the media is piecewise constant to see that for this case we can achieve considerably better results. In [4] pth order accurate SBP op-erators D2 for the second derivative was constructed. Such an operator has

the properties

D2 = H−1(−M + BS), H = HT > 0 (18)

M = MT > 0, B = diag(−1, 0, . . . , 1)

and S approximates the first derivative at the boundaries. Further D2 was

constructed to have the narrowest interior stencil possible. When H is diag-onal D2 is called a pth order accurate narrow - diagonal SBP operator for

the second derivative. These operators has the properties that in contrast to using D1D1 as an approximation to the second derivative the highest

fre-quency that can reside on the grid, the π - mode, is damped by D2 whereas

it is not by D1D1. Also, the interior stencil is wider in D1D1, see Figure (7).

0 5 10 15 20 0 2 4 6 8 10 12 14 16 18 20 22 nz = 243 Non − zero elements of D₁ D₁

(a) Non - zero elements of D1D1

0 5 10 15 20 0 2 4 6 8 10 12 14 16 18 20 22 nz = 147 Non − zero elements of D₂

(b) Non - zero elements of D2

Figure 7: The stencil width is narrower in the narrow diagonal SBP operator

D2 compared to D1D1. In this figure 6th order SBP operators are depicted.

As before we consider a domain consisting of two sub - domains coupled at the east respective west boundaries. We denote by v the solution on the left block and by w the solution on the right block. Let

I1 = vE− wW = 0, I2 = (b1(BS)xv)E + (b2(BS)xw)W = 0.

I3 = vtE− wtW = 0, I4 = vttE − wttW = 0.

The semi - discretization with homogeneous Neumann boundary conditions looks like vtt = b1D2xv + b1D2yv wtt = b2D2xw + b2D2yw + τ H_x−1iEI1 − τHx−1iWI1 +1 2H −1 x b1((BS)x)TiEI1 − 1 2H −1 x b2((BS)x)TiWI1 −1 2H −1 x iEI2 − 1 2H −1 x iWI2 (19) + σH_x−1iEJ3 − σHx−1iWI3 + ρH_x−1iEJ4 − ρHx−1iWI4 − H−1 x iWb1(BS)xv − Hx−1iEb2(BS)xw − H−1 y iNb1(BS)yv − Hy−1iNb2(BS)yw − H−1 y iSb1(BS)yv − Hy−1iSb2(BS)yw

The proof of the stability of the above scheme requires the following corollary to lemma 2.3 of [4]

Corollary 2 For D2 and u a 2 - dimensional grid function uT(bM ⊗ H)u = hα b((BS)xu) T WH((BS)xu)W + h α b((BS)xu) T EH((BS)xu)E + uT( ˜bM ⊗ H)u

where ˜M is symmetric and positive definite and α a positive constant inde-pendent of h.

Proof It is proved in [5] that vTbM v = hα_b(bBSv)2₀+ hα_b(bBSv)2_N+ vTbM v˜

holds for any 1 - dimensional grid function v. Let now u(i) be the vector of

values of u at yi in the x - direction. We get

uT(bM ⊗ H)u = ΣiHiiuT(i)bM u(i)

= hα

bΣiHii

(

Second - order Fourth - order Sixth - order

0.4 0.2508560249 0.187815026

Table 4: The values of α in lemma 2 for different operators

The constant α of lemma 2 is dependent on the order of accuracy of the second derivative SBP operator used, values for α was derived in [5] and are listed in Table 4. Now, we get as corollary to lemma 2.4 in [5].

Corollary 3 If τ ≤ −b1+b2

4αh , σ ≤ 0 and ρ ≤ 0, then the scheme (19) yields a

strictly stable method for solving the model problem in Cartesian coordinates.

Proof:

Follow the reasoning of [5] analogous to the 2 dimensional case and use corollary 2 in the place of lemma 2.3 of [5] △

We introduce the penalty parameter ΓCa by

τ = ΓCa(−

b1+ b2

4αh )

We can now perform the same test of accuracy and rate of convergence as in the first case of section 4 above. Table 5 should be compared with table 1 and Table 6 should be compared with Table 2. The improvement in using the narrow operator D2instead of D1D1is striking. Now, by these comparisons if

4th order interior stencil

ΓCa= 1 ΓCa= 1.2 ΓCa= 5

N log₁₀er(N ) q(N) log10er(N ) q(N) log10er(N ) q(N)

26 -1.81 - -2.10 - -2.47

-51 -2.45 2.20 -3.49 4.78 -3.68 4.14

101 -3.18 2.45 -4.90 4.73 -4.95 4.29

201 -3.93 2.50 -6.13 4.14 -6.20 4.19

6th order interior stencil

ΓCa= 1 ΓCa= 1.2 ΓCa= 5

N log₁₀er(N ) q(N) log10er(N ) q(N) log10er(N ) q(N)

26 -1.73 - -2.17 - -2.51

-51 -2.96 4.17 -4.09 6.52 -4.53 6.88

101 -4.28 4.47 -6.07 6.70 -6.58 6.94

201 -5.63 4.51 -8.15 6.99 -8.45 6.23

Table 5: Rate of convergence and accuracy for the scheme (19) for different values of ΓCa without the extra conditions.

ΓCa = 8, σ =−4, ρ = −4

4th order interior stencil 6th order interior stencil N log₁₀er(N ) q(N) log10er(N ) q(N)

26 -2.47 - -2.52

-51 -3.68 4.13 -4.51 6.82

101 -4.95 4.28 -6.60 7.05

201 -6.20 4.18 -8.47 6.26

Table 6: Rate of convergence and accuracy for the scheme (19) with the extra conditions.

References

[1] H - O Kreiss, J Oliger, Comparison of accurate methods for the

integra-tion of hyperbolic equaintegra-tions, Tellus, No 24 (1972)

[2] B Gustafsson, H - O Kreiss, J Oliger, Time dependent problems and

diifference methods, Wiley 1995

[3] B Strand, Summation by parts for finite difference approximations for

d/dx, J. Comp Phys 110 (1994)

[4] K Mattsson, J Nordström, Summation by parts operators for finite

[5] K Mattsson, F Ham, G Iaccarino, Stable and accurate wave - propagation

in discontinuous media, J. Comp Phys 227 (2008)

[6] L N Trefethen, D Bau, III Numerical Linear algebra, SIAM (1997) [7] K Mattsson, M Carpenter, Stable and accurate interpolation operators

for high - order multi - block finite difference methods, SIAM Journal on

Scientific Computing, Vol 32, No 4 (2010)

[8] K Mattsson, Boundary procedures for summation - by - parts operators, Journal of scientific computing Vol 18, No 1 (2003)

[9] K Mattsson, J Nordström, High order finite difference methods for wave

propagation in discontinuous media, J. Comp Phys 220 (2006)

[10] D Appelö, N A Petersson, A stable finite difference method for the

elas-tic wave equation on complex geometries with free surfaces,

(a) t = 0

(b) t = 0.25