Stable and high order accurate difference methods for the elastic wave ...

Journal of Computational Physics 279 (2014) 37–62

Contents lists available at ScienceDirect

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

Stable and high order accurate difference methods for the elastic wave equation in discontinuous media Kenneth Duru a , Kristoffer Virta b a b

Department of Geophysics, Stanford University, Stanford, CA, United States Division of Scientific Computing, Uppsala University, Uppsala, Sweden

a r t i c l e

i n f o

Article history: Received 5 December 2013 Received in revised form 26 August 2014 Accepted 28 August 2014 Available online 6 September 2014 Keywords: Elastic waves Interface waves Reflected waves Refracted waves High order accuracy Wellposedness Time-stability SBP–SAT

a b s t r a c t In this paper, we develop a stable and systematic procedure for numerical treatment of elastic waves in discontinuous and layered media. We consider both planar and curved interfaces where media parameters are allowed to be discontinuous. The key feature is the highly accurate and provably stable treatment of interfaces where media discontinuities arise. We discretize in space using high order accurate finite difference schemes that satisfy the summation by parts rule. Conditions at layer interfaces are imposed weakly using penalties. By deriving lower bounds of the penalty strength and constructing discrete energy estimates we prove time stability. We present numerical experiments in two space dimensions to illustrate the usefulness of the proposed method for simulations involving typical interface phenomena in elastic materials. The numerical experiments verify high order accuracy and time stability. © 2014 Elsevier Inc. All rights reserved.

1. Introduction The elastic wave equation governs the propagation of seismic waves resulting from earthquakes and other seismic events. Other applications include waves in plates, beams, solid structures and crystals. In the general setting the material parameters describing the media are piecewise smooth functions with discontinuities at material boundaries. This is especially true in seismological problems due to the layered structure of the Earth. Large contrasts in media parameters can also be found in solid mechanics devices composed of different materials in welded contact. The presence of material discontinuities gives rise to reflection and refraction phenomena. For example earlier than expected arrivals of seismic waves, so called refraction arrivals, have been shown to exist in the presence of a material discontinuity [9,24]. Similarly to the Rayleigh surface wave [25] an analogous Stoneley interface wave may under certain circumstances travel along media interfaces [27]. Also, as in the case of the traction free boundary a distinguishing characteristic of wave–interface interaction is that mode conversion occurs. That is, an incident wave, either pressure or shear, is converted into two wave types, pressure and shear on reflection and refraction. More details on the theory of interaction of wave modes generated by material boundaries or discontinuous media can be found in [1,11]. In this paper, we consider the elastic wave equation in a bounded and discontinuous medium. We are particularly interested in accurate and efficient numerical simulation of wave phenomena arising from the interfaces of internal material boundaries. High order accurate schemes have been demonstrated to be more efficient than lower (first and second) order accurate methods for wave propagation problems in smooth domains. However, the presence of boundary conditions and discontinuities in media parameters makes the design of stable and accurate numerical methods challenging. There are two main approaches for discretizing the elastic wave equation. A common approach is to rewrite the elastic wave equation http://dx.doi.org/10.1016/j.jcp.2014.08.046 0021-9991/© 2014 Elsevier Inc. All rights reserved.

38

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

as a first order system (using for example the velocity-stress formulation) prior to numerical approximations. This is in part due to maturity of the theory and numerical methods for first order hyperbolic systems [13,26,29]. Another approach discretizes directly the elastic wave equation in a second order form where the unknown variables are the displacement fields [7,23,36–38]. In the finite difference community, staggered grid approximations [32–35] for first order hyperbolic system are quite popular because of their low dispersion errors. However, the introduction of boundary conditions for high order accurate approximations are cumbersome using the staggered grid approach. It is important to note that narrow stencil finite difference operators for second derivatives [7,17,23] also have optimal dispersion properties, equivalent to staggered grid methods. A recent and promising strategy for numerical treatments of second order systems combining narrow stencil summationby-parts (SBP) finite difference operators with the use of ghost points to enforce boundary and interface conditions are described in [16,23,36–38]. In [16] analysis of numerical errors for Rayleigh surface waves in almost incompressible materials was presented. The result is that, to ensure efficient numerical treatments higher order accurate schemes (at least 4th order accuracy) are essential. In [23] the focus was put on the stable treatment of hanging nodes in a grid refinement and accurate discretization of singular source terms. Numerical experiments in a discontinuous medium was also presented with the results showing good agreement with a semi-analytical solution. In a recent work [7], the elastic wave equation in second order form with smooth material coefficients and subject to general (including traction free and Dirichlet) boundary conditions was considered. For the first time, a sixth order accurate and time-stable approximation applicable to the elastic wave equation in heterogeneous and complex geometry was derived by using SBP finite difference operators to approximate spatial derivatives and weak enforcements of boundary conditions. The results of [7] and the analysis in [30] demonstrate the efficiency of our high order methods in accurate simulation of surface waves in complex elastic media. The present study should be seen as a direct continuation of the work in [7]. We consider the elastic wave equation in second order form, and seek an approximation of the displacement field due to an initial disturbance as a function of time. A possible approach to ensure high order accuracy across discontinuous interfaces is to decompose the problem into multiple subdomains, where each block is characterized by smooth coefficients. Reflection and transmission (interface) conditions are imposed at the interfaces. Then each subdomain can be discretized using a high order accurate scheme. The local subdomains are patched together into the global domain in a stable manner. This approach is often referred to as domain decomposition. In the finite difference setting, a stable and accurate domain decomposition strategy can be achieved using the multi-block structured approach [3]. The multi-block structured scheme suggests: decompose the domain into a number of smooth and structured meshes in a multi-block fashion; discretize the sub-domains using high order accurate SBP operators; and impose the interface conditions weakly using a penalty technique referred to as simultaneous approximation term (SAT) [2]. The critical aspect of this procedure is choosing penalty parameters and ensuring numerical stability. For first order systems, the SBP–SAT approach for interface treatments is rather well understood [3,13]. It has also been recently extended to the scalar wave equation in second order form [6,21] on Cartesian grids. In [17] SBP operators for second derivatives with variable coefficients are derived. These operators are used in [31] to construct stable and high order accurate numerical approximations for the scalar wave equation in complex geometries. The main challenge for the elastic wave equation is that the spatial energy operator has a nontrivial null space. Therefore, the technique developed in [31] for the scalar wave equation cannot be easily extended to the elastic wave equation. In addition, the elastic wave equation in second order form contains second derivative terms such as ∂/∂ x(a∂ u /∂ x) and mixed derivative terms like ∂/∂ x(b∂ u /∂ y ), where a > 0, b ∈ R. The second and mixed derivative terms appear in the interface conditions as a linear combination of normal and tangential derivatives on the boundaries. In a finite difference setting this poses a numerical challenge. To discretize these terms, approximations of both second and mixed derivatives terms, as well as approximations of normal and tangential derivative terms on the boundaries are required. Then certain compatibility conditions of the approximations need to be fulfilled. We use SBP approximations of first and second derivatives that are either compatible or fully compatible (the definitions will become clear in Section 3). The interface conditions are imposed weakly using penalties. By deriving penalty parameters and corresponding discrete energy estimates we prove time stability. By using the time-stepping scheme derived in [8] we obtain a fully discrete high order accurate and time-stable scheme. The main result of the paper is the theoretical high order accuracy of our method across discontinuous interfaces and the ‘bulletproof’ of numerical stability. Using an unstructured mesh greatly simplifies the handling of complicated geometries, in particular if hanging nodes are allowed. Methods with this approach include discontinuous Galerkin [14] and finite element [4] discretizations of the elastic wave equation. However, generating a high quality unstructured mesh can however be a non-trivial problem in itself. Unstructured discretizations also require extra bookkeeping and additional memory to keep track of the connectivity of the grid. In this paper we consider problems where multi-block curvilinear grids can be generated efficiently. We believe that the schemes described in this paper provide the opportunity for the development of powerful (fast and accurate) solvers for various wave propagation problems encountered in geophysics and solid mechanics problems. The rest of the paper is outlined as follows. In Section 2 the equations and interface conditions are introduced in the general setting. Section 3 introduces the necessary definitions used to describe the discretization. Stability of the numerical scheme is proved in Section 4. In Section 5 numerical experiments are presented. The purpose of the experiments are twofold: One, to verify stability and accuracy, two, to illustrate the usefulness of the proposed method on simulation of phenomena typical to wave motion in discontinuous elastic media. In particular, numerical examples of mode conversion and the Stoneley interface wave are given. Section 6 concludes and mentions future work.

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

39

2. Equations of linear elasticity In this section, we present the equations of linear elastodynamics in discontinuous media. First we consider Cartesian coordinates. Second we consider curvilinear coordinates. 2.1. Elastic waves in rectangular Cartesian coordinates Let the Cartesian coordinates be described by (x, y ) with x denoting the horizontal direction and y denoting the vertical direction. Consider elastic waves propagating in two (upper and lower) heterogeneous anisotropic elastic half-planes, whose principal axes coincide with (x, y ). The half-planes are in a welded contact at y = 0. That is, slip is not allowed at the interface. The displacement fields u = (u 1 , u 2 ) T in the lower half-plane y < 0 and v = ( v 1 , v 2 ) T in the upper half-plane y > 0 are governed by:

        ∂ ∂ ∂ ∂ ∂ 2u ∂u ∂u ∂u T ∂u + + + , y < 0, = A B C C ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x ∂t2         ∂ ∂ ∂ ∂ ∂ 2v ∂v ∂v ∂v ∂v + + + , y > 0. ρ 2 = A B C C T ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x ∂t

ρ

(1) (2)

Here, the densities ρ > 0, ρ  > 0 and the coefficients matrices A, B, C and A  , B  , C  define the media. The systems (1), (2) are symmetric hyperbolic. Thus the coefficient matrices are real and satisfy

  (kx , k y ) = k2x A + k2y B + kx k y C + C T , 

T

 =  ≥ 0,

 =

T

≥ 0,

    (kx , k y ) = k2x A  + k2y B  + kxk y C  + C  T , 2

∀(kx , k y ) ∈ R .

Note that since the matrices  ,   are symmetric positive semi-definite the eigenvectors of

(3) (4)

 ,   must be orthogonal and

the eigenvalues real and nonnegative. As a matter of fact, there must be at least one nontrivial eigen-pair corresponding to each of the matrices  ,   . At y = 0, the traction vectors and displacement vectors are continuous,

 



T n = T n,

n=

0 1

u = v.

,

(5)

The stress tensors T , T  are given by

  ∂u ∂u ∂u ∂u , T = A +C ,B + CT ∂x ∂y ∂y ∂x

  ∂v ∂v ∂v ∂v T  = A . + C  , B + C T ∂x ∂y ∂y ∂x

(6)

Thus, by using (5) and (6) we have

B

∂u ∂u ∂v ∂v + CT = B + C T , ∂y ∂x ∂y ∂x

u = v,

at y = 0.

(7)

Let P , P  denote the 4 × 4 potential energy matrices defined by



P=

A CT

C B



,

P =



A

C

C T

B



(8)

.

The potential energy matrices P and P  are symmetric and positive semi-definite matrices, see [7]. We introduce the standard L 2 -scalar product with the corresponding norm





vT udx,

(u, v) =

u2 =

Ω

uT udx.

(9)

Ω

The elastic energy in each medium can be defined by [7]

    ∂ u T  ∂ u   √ ∂ u 2 ∂x ∂x ρ  + P ∂ u dxdy , ∂u ∂t ∂y ∂y

Eu (t ) :=  

Ω1

     ∂ v 2   + ρ E v (t ) :=   ∂t  Ω2

∂v ∂x ∂v ∂y

T

P



∂v ∂x ∂v ∂y

(10)

 dxdy .

(11)

The total energy is given by the sum of the energies

E(t ) = Eu (t ) + E v (t ).

(12)

40

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Since we consider half-planes there are no boundaries in the x-direction, and at | y | → ∞. It is then straightforward to show that the elastic wave equations (1), (2) with the interface condition (7) satisfy

E(t ) = E(0),

∀t > 0.

(13)

The elastic energy (12) is conserved. Thus, the systems (1) and (2) with the interface conditions (7) is well-posed. In an orthotropic anisotropic media whose principal axis coincides with the x, y axis, the coefficient matrices have a simple structure of the form



A=

c 11 0



A =

0 c 33





B=

,  0 , 

 c 11

0

B =

c 33

c 33 0



0 c 22

 c 33

0





C=

, 

0

C =

,

c

0 c 33



c 12 0

0 c

22

33



 c 12

(14)

, 

0

(15)

.

 , c  , c  , c  are the elastic coefficients Here, c 11 , c 12 , c 22 , c 33 are the elastic coefficients in the upper half-plane and c 11 12 22 33     can be negative (for certain in the lower half-plane. The coefficients c 11 , c 22 , c 33 , c 11 , c 22 , c 33 are positive and c 12 , c 12 materials). For orthotropic media the elastic coefficients satisfy

c 33 ≥ 0,

c 11 , c 22 > 0,

  c 11 , c 22 > 0,

2 c 11 c 22 − c 12 > 0,

By using the relation (16) we also have



 c 33 ≥ 0,

  2 c 11 c 22 − c 12 > 0.

(16)



2 2 < (c 11 + c 22 )2 , 0 ≤ (c 11 − c 22 )2 + 4c 12 = (c 11 + c 22 )2 − 4 c 11 c 22 − c 12



  0 ≤ c 11 − c 22

2

          2  2 2  2 < c 11 + c 22 + 4c 12 = c 11 + c 22 − 4 c 11 c 22 − c 12 .

 = c  = 2μ + λ , c  = μ , c  = λ , and In isotropic elastic media we have c 11 = c 22 = 2μ + λ, c 33 = μ, c 12 = λ, c 11 22 33 12 −μ < λ < ∞, −μ < λ < ∞. In the special case of a solid/fluid interface we must have either μ = 0 or μ = 0.

2.2. Elastic waves in curvilinear coordinates Often the underlying media can be much more complicated such that material interfaces are no longer simple planar boundaries. Efficient numerical treatments of curved interfaces require numerical approximations on curvilinear grids. By splitting the domain and introducing the smooth grid transformation (x(q, r ), y (q, r )) ↔ (q(x, y ), r (x, y )), for each block, the elastic wave equation can be written in curvilinear coordinates

        ∂ ∂u ∂ ∂u ∂ T ∂ u ∂ ∂u ∂ 2u + + + , r < r0 , = A B C C ∂q ∂q ∂r ∂r ∂q ∂r ∂r ∂q ∂t2         ∂  ∂ v ∂  ∂ v ∂  T ∂ v ∂  ∂ v ∂ 2v + + + , r > r0 , ρ  2 = A B C C ∂q ∂q ∂r ∂r ∂q ∂r ∂r ∂q ∂t

ρ

where

and

(17) (18)

    A = J q2x A + q2y B + q x q y C + q x q y C T , B = J r x2 A + r 2y B + r x r y C + r x r y C T ,   C = J qxrx A + q y r y B + qxr y C + q y rx C T , ρ = J ρ

(19)

    A  = J q2x A  + q2y B  + q x q y C  + q x q y C  T , B  = J r x2 A  + r 2y B  + r x r y C  + r x r y C  T ,   C  = J qx rx A  + q y r y B  + qxr y C  + q y rx C  T , ρ  = J ρ  .

(20)

The Jacobian of the transformation is

J = xq y r − xr yq > 0, and the metric relations are given by

J qx = yr ,

J q y = −xr ,

J rx = − yq ,

J r y = xq .

Note that since the elastic wave equations (1), (2) are symmetric hyperbolic the transformed coefficient matrices A, A  and

B, B  are symmetric positive semi-definite. Thus we have



A=

S TA

B=

S TB



(a)

θp

0 (b)

θp

0

0



(a)

θs

0 (b)

θs



S A,



A =

S TA

SB,



S TB



A =

(a )

θp

0 (b )

θp

0

0



(a )

θs

0 (b )

θs





(b )

(b )

S A,

θ p(a) , θs(a) , θ p(a ) , θs(a ) ≥ 0.

SB,

θ p , θs , θ p , θs

(21)

 (b)

(b)

≥ 0.

(22)

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

41

In orthotropic elastic media the eigenvalues are (a)

θk =

(a )

θk

θk

 (c 11 + c 33 )q2x + (c 22 + c 33 )q2y 2      J  2 (c 11 + c 33 )q2x + (c 22 + c 33 )q2y − 4 c 11 c 33 q4x + c 22 c 33 q4y + c 11 c 22 + c 33 ± − (c 12 + c 33 )2 q2x q2y ,

2  2     J    = c 11 + c 33 q x + c 22 + c 33 q2y 2            J    q2 + c  + c  q2 2 − 4 c  c  q4 + c  c  q4 + c  c  + c  2 − c  + c  2 q2 q2 ± c 11 + c 33 x x y 33 22 33 y 11 33 x 22 33 y 11 22 12 33 2 (23)

θk(b) =

(b )

J

2

J

 (c 11 + c 33 )r x2 + (c 22 + c 33 )r 2y 2 2      J  2 (c 11 + c 33 )r x2 + (c 22 + c 33 )r 2y − 4 c 11 c 33 r x4 + c 22 c 33 r 4y + c 11 c 22 + c 33 ± − (c 12 + c 33 )2 r x2 r 2y ,

2 2    2 J     = c 11 + c 33 r x + c 22 + c 33 ry 2            J    r2 + c + c r2 2 − 4 c c r4 + c c r4 + c c + c 2 − c + c 2 r2r2 , ± c 11 + c 33 x x y 33 22 33 y 11 33 x 22 33 y 11 22 12 33 2 (24)

where k = p , s. The eigenvalues with positive square roots (k = p) correspond to the quasi-P waves and the eigenvalues with negative square roots (k = s) correspond to the quasi-S waves. In linear isotropic elastic media, the eigenvalues are

  

 (a) θk = μ J q2x + q2y , (λ + 2μ) J q2x + q2y ,      

(a ) θk = μ J q2x + q2y , λ + 2μ J q2x + q2y ,   

 (b) θk = μ J r x2 + r 2y , (λ + 2μ) J r x2 + r 2y ,      

 θk(b ) = μ J r x2 + r 2y , λ + 2μ J r x2 + r 2y . Similarly, the eigenvalues above correspond to the P- and S-waves propagating in the media. The new coordinates q, r are regular Cartesian grids with a planar interface at r = 0. At the interface at r = 0, traction vectors and displacement vectors are continuous,

T n = T  n,

1 n= 2 r x + r 2y





rx

u = v.

,

ry

(25)

Thus we have

    ∂ u T ∂ u 1 ∂ v  T ∂ v = , B +C B +C ∂r ∂q ∂r ∂q r2 + r2 r2 + r2 1



x

y

x

u = v.

(26)

y

Similarly, in the transformed coordinates we introduce the 4 × 4 symmetric semi-positive definite [7], potential energy, matrices P, P  defined by

P=

 A CT

 C , B

P =

  A C T

C B



where

P = S P ST ,

P = S PST ,

(27)

 S=

qx

rx

qy

ry



 ⊗

1 0 0 1

 ,

det ( S ) =

1 J2

and ⊗ is the Kronecker product defined in Appendix A. Notice also from (28) we have

P = P T ≥ 0 ⇐⇒ P = P T ≥ 0,

P  = P  T ≥ 0 ⇐⇒ P = P  T ≥ 0.

= 0,

(28)

42

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

The elastic energy in each medium can be defined by [7]

    ∂ u T  ∂ u   ∂ u 2 ∂q ∂q   ˜ Eu (t ) :=  ρ + P dqdr , ∂u ∂u ∂t  1 Ω

∂r

2 Ω

∂r

(29)

∂r

    ∂ v T  ∂ v   ∂ v 2 ∂q ∂q   + E v (t ) :=  ρ P dqdr .  ∂v ∂v ∂t 

(30)

∂r

The total energy is given by the sum of the energies

E(t ) = Eu (t ) + E v (t ).

(31)

Analogously, by ignoring boundaries in the q-direction and at |r | → 1, it is straightforward to show that the elastic wave equations (17), (18) in curvilinear coordinates with the interface condition (7) satisfy

E(t ) = E(0),

∀t > 0.

(32)

The sum of the elastic energy (31) is conserved. 3. Spatial discrete approximations In this section, we introduce the spatial operators used in this study. We begin with the approximations of the first and second spatial derivatives of a scalar variable in one space dimension, and extend the spatial operators to vector variables and two space dimensions. We end the section by deriving a stable numerical interface treatment for elastic wave equations. A numerical method for the outer boundaries was recently considered in [7]. 3.1. One space dimension Let the unit interval r ∈ [0, 1] be discretized with the uniform grid

r j = ( j − 1)hr ,

hr =

1 Nr − 1

j = 1, · · · , N r .

,

The grid function at r j is denoted v j (t ) and the semi-discrete solution vector is v(t ) = [ v 1 (t ), v 2 (t ), v 3 (t ), . . . , v N (t )] T . The finite difference matrices are denoted

Dr ≈

∂ , ∂r

(a)

D rr ≈

  ∂ ∂ . a(r ) ∂r ∂r (a)

Definition 1. The matrices D r , D rr are (one dimensional) summation-by-parts (SBP) operators [17] if they satisfy the following properties

D r = H r−1 Q r ,

 −1

(a)

D rr = H r

Q rT + Q r = E Rr − E Lr ,



(a)

− M r + ( E Rr − E Lr )aS r ,

(a)

v T M r v ≥ 0,

H r = H rT ,

(a)

Mr



(a)  T

= Mr

(33)

,

v T H r v > 0.

(34)

Here, a = diag(a(r1 ), a(r2 ), · · · , a(r Nr )) and E Lr = diag(1, 0, 0, . . . , 0), E Rr = diag(0, 0, 0, . . . , 1) pick out the left and right (r ) boundary terms. The matrix H r is diagonal with H r j j = h j j > 0, thus defines a discrete norm. The matrix Q r is almost (a)

skew-symmetric, M r

(a)

is called the symmetric part of D y y and S r is a consistent approximation of the first derivative ∂/∂ r (a)

on the boundaries. The interior scheme of D rr is a centered narrow stencil approximation, as opposed to applying the first derivative operator D r twice. (a)

Definition 2. Let D r and D rr , defined in (33), (34) above, denote SBP operators approximating ∂/∂ r, ∂/∂ r (a(r )∂/∂ r ) with 2s interior order of accuracy and s accurate boundary closures. Let the corresponding diagonal norm be denoted by H r . The (a) operators D r and D rr are compatible SBP operators if (a)

Mr

(a)

= D rT H r aD r + R r ,

(a)

Rr

 (a) T = Rr ,

(a)

v T R r v ≥ 0.

(35)

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

43

We can write, see [17], (a)

Rr

=

s 

T

α (k) D (k) C ((ak)) D (k) .

(36)

k =1

(k)

Here, α (k) > 0 are weights, D (k) are difference operators and C (a) are diagonal matrices with nonnegative entries. For constant coefficients, narrow stencil compatible SBP operators can be found in [20–22]. The compatible SBP operators for variable coefficients are recently derived in [17]. (a) Let D rr denote an SBP operator [17,22] approximating ∂/∂ r (a(r )∂/∂ r ) with 2s interior order of accuracy and correspond(a) (a) ing to a diagonal norm defined by H r . An important term in the derivation of truncation errors is ( H r−1 M r U ) where M r (a)

is the symmetric part of D rr . By Taylor expansions we obtain

     −1 (a)  du d a(r1 ) du (r1 ) + (r ) − H r Mr U 1 = a (r1 ) + O hrs , dr



dr

h11 dr



    du d − H −1 M r(a) U j = (r j ) + O hrs , 2 ≤ j ≤ N s , a dr dr      du d (a)  (r j ) + O hr2s , N s + 1 ≤ j ≤ N r − N s , a − H −1 M r U j = dr dr      du d (a)  (r j ) + O hrs , N r − N s < j ≤ N r − 1, a − H −1 M r U j = dr dr      −1 (a)  du d a(r N ) du (r Nr ) − (r ) r a (r Nr ) + O hrs . − H Mr U N = dr

r

dr

h N r N r dr

(37)

Here, N s denotes the number of grid points where lower order accurate finite difference stencils are used. Apparently, in (a) order for D rr to yield sth order accuracy at the boundary grid points, the accuracy of the first derivative operator S r must be at least (s + 1). Note that the first derivative operator S r is not unique. We can therefore replace S r by D r in (34) without a significant effect on accuracy. (a)

Definition 3. Consider the SBP operators D rr and D r defined in (34), (33). The operators are called fully compatible SBP operators if for the norm H r , (35) holds and S r = D r . (a)

Lemma 1. Let D rr ≈ ∂/∂ r (a(r )∂/∂ r ) denote the diagonal norm SBP operator, defined in (34), with 2s interior order of accuracy and s order of accuracy close to the boundaries. If D r ≈ ∂/∂ r is a diagonal norm SBP operator of the same order of accuracy, then the operator



¯ rr = H r−1 − M r D (a)

(a)

 + ( E Rr − E Lr )aD r ,

(38)

is a diagonal norm SBP operator with 2s interior order of accuracy and s order of accuracy close to the boundaries, with the exception of the boundary points, j = 1, N r , where the order of accuracy is s − 1. Proof. Introduce the smooth function u (r ) and denote U (r j ) the restriction of the function on the grid with N r number of

¯ rr U ) j at the grid point r j yields grid points. Taylor expansions of the expression ( D (a)



¯ rr U D (a)



     p   −1  du d (a) − (r j ) + O hr j . = H M + ( E − E ) aD U = a Rr Lr r r r j j dr

dr

(39)

We will use (37) and the Taylor expansion of the boundary stencil for the first derivative operator ( D r U ) j . Note that at the boundary points the accuracy of the difference operator ( D r U ) j is s. Thus, we have

   (r1 ) + O hrs−1 , dr dr    (a)    ¯D rr U = d a du (r j ) + O hrs , 2 ≤ j ≤ N s , j dr dr    (a)    ¯ rr U = d a du (r j ) + O hr2s , N s < j ≤ N r − N s , D j dr dr    (a)    du d ¯ rr U = (r j ) + O hrs , N r − N s < j ≤ N r − 1, D a j dr dr    (a)    du d ¯ rr U (r Nr ) + O hrs−1 . D = a N 

¯ rr U D (a)



= 1

r

d



dr

a

du

dr

This completes the proof of the lemma.

2

(40)

44

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

For constant coefficients, fully compatible SBP operators can be found in [19]. To the best of our knowledge, fully compatible SBP operators for second derivatives with variable coefficients are not yet available. However, by Lemma 1 above, fully compatible SBP operators for variable coefficients can be adapted from the compatible SBP operators derived in [17], by using S r = D r in (34). We have used these adapted fully compatible SBP operators, for second derivatives with variable coefficients, in the numerical experiments performed in Section 5.4 of this paper. 3.2. Two space dimensions Consider the transformed elastic wave equation in a two layered elastic medium Ω = Ω1 ∪ Ω2 , where Ω1 = [0, 1] × [0, 1] and Ω2 = [0, 1] × [−1, 0]. Note that on rectangular Cartesian grids x = q, y = r, and the metric coefficients and Jacobian are q x = 1, r x = 0, q y = 0, r y = 1 and J = 1. Here, we will focus on the numerical treatment of the interface conditions (7) at r = 0. A numerical method for the outer boundaries was recently considered in [7]. To begin with, in the upper and lower layers we introduce a uniform grid defined by constant spatial steps hq = 1/( N q − 1), hr = 1/( N r − 1). Here, N q is the number of grid-points used in the q-direction and N r is the number of grid-points used in the r-direction. We denote the grid function [ui , j ]. The grid function [ui , j ] is stacked as vector of length N q N r



u = (u )11 , (u )12 , . . . , (u ) N q N r

T

.

The semi-discrete solution is a vector of length 2N q N r . All scalar two space dimensional material parameters a(x, y ) are rearranged as diagonal matrices as follows





a = diag a(q1 , r1 ), a(q1 , r2 ), . . . , a(q Nq , r Nr ) . In order to enable compact notation we will use the Kronecker product defined in Appendix A. Let, a11 , a12 , a21 , a22 denote diagonal matrices of size N q N r × N q N r for the material parameters. Correspondingly, A, B, C, C T are block matrices on the form

A=



a11 a12 a12 a22



B=

,



b11 b12



b12 b22

,

C=

  c11 c12 , c21 c22

CT =

  c11 c21 . c12 c22

Discrete spatial operators in two space dimensions are defined by

Dq = I 2 ⊗ D q ⊗ I r ,



( A)

Dqq =

( a11 )

( a12 )

Dqq

Dqq

( a12 )

Dr = I 2 ⊗ I q ⊗ D r ,



Dqq

( B)

Drr =

,

( a22 )

Dqq



( b11 )

Sq = I 2 ⊗ S q ⊗ I r ,

 ( b12 )

Drr

Drr

Drr

Drr 22

( b12 )

( b )

Sr = I 2 ⊗ I q ⊗ S r ,

(41) (42)

.

Note that



   + I 2 ⊗ ( E Rq − E Lq ) ⊗ I r A( I 2 ⊗ S q ⊗ I r ) ,     ( A) ( A) Drr = ( I 2 ⊗ I q ⊗ H r )−1 −Mr + I 2 ⊗ I q ⊗ ( E Rr − E Lr ) A( I 2 ⊗ I q ⊗ S r ) , ( A)

( A)

Dqq = ( I 2 ⊗ H q ⊗ I r )−1 −Mq

(43)

where

  ( A) = I 2 ⊗ D qT ⊗ I r ( I 2 ⊗ H q ⊗ I r ) A( I 2 ⊗ D q ⊗ I r ) + Rq ,   ( B) ( B) Mr = I 2 ⊗ I q ⊗ D rT ( I 2 ⊗ I q ⊗ H r ) B( I 2 ⊗ I q ⊗ D r ) + Rr , ( A)

Mq

 T = Rq(A) ≥ 0,  ( ( B) B)  T Rr = Rr ≥ 0. ( A)

Rq

Here, S q , D q , E Rq , E Lq are square matrices with the dimension N q × N q and S r , D r , E Rr , E Lr are square matrices with the dimensions N r × N r . The matrices S q , S r represent one-sided difference approximations of the normal derivatives, defined in the SBP operator (34) for second derivatives. The matrices D q , D r are (one dimensional) SBP operators (33) approximating the first derivatives. The matrices E Rq = diag(0, . . . , 1), E Lq = diag(1, 0, . . . , 0) and E Rr = diag(0, . . . , 1), E Lr = diag(1, 0, . . . , 0) pick out the boundary terms. By using (36) we can rewrite the remainder operators as ( A)

Rq

=

r  j =1

( B)

Rr

=

r  j =1

( j) (A)

( j) (B)



T







T





α ( j) I 2 ⊗ D ( j) ⊗ I r C(( Aj)) I 2 ⊗ D ( j) ⊗ I r , 

α ( j) I 2 ⊗ I q ⊗ D ( j) C(( Bj)) I 2 ⊗ I q ⊗ D ( j) ,

where C , C are nonnegative block diagonal matrices.

(44)

(45)

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

45

3.3. Numerical interface treatment We will now define the discrete spatial approximations for the right hand sides of (1) and (2) ( A) ( B) P = Dqq + Drr + Dq CDr + Dr C T Dq , ( A )

(46)

( B )

P  = Dqq + Drr + Dq C Dr + Dr C T Dq .

(47)

To extract the stresses on the boundaries we will use



q2x + q2y Tq = ASq + CDr ,



q2x + q2y Tq = A Sq + C Dr ,



r2x + r2y Tr = BSr + C T Dq ,



r2x + r2y Tr = B Sr + C T Dq .

Note also that the restriction of the discrete operators Tq on the vertical boundaries, and Tr on the horizontal boundaries applied to the grid function yield approximations of the normal stresses on the boundaries. In order to enforce the interface conditions (7), we introduce two interface operators BI ,  B I defined by

⎧    ⎨0 u BI = ui , N y − vi ,1 v ⎩ ij vi ,1 − ui , N y ⎧    ⎨0 u  BI = ui , N y − vi ,1 v ⎩ ij ui , N y − vi ,1

if j = { N y , N y + 1}, if j = N y ,

(48)

if j = N y + 1, if j = { N y , N y + 1}, if j = N y ,

(49)

if j = N y + 1.

The interface operators compute the difference of the solutions on the interface and projects them to the corresponding equations. At two grid-points collocated on the interface, the interface operators BI,  B I have a compact matrix notation,

⎡

0 ⎢0

⎢ ⎢ ⎢ ⎢ ⎢0 ⎢ ⎢0 BI = ⎢ ⎢0 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎣

0 0

0 0

..

..

. ... ... ... ...

.

0 0 0 0

..

.

0 0 1 −1 −1 1 0 0

0 0 0 0

... ... ... ... .. .. . . 0 0

0 0

DCq =

,

 T

C Dq 0

Hqr = I 2 ⊗ H q ⊗ H r , 0

C T Dq

⎛

 ,

Tr = ⎝

⎢ ⎢ ⎢ ⎢ ⎢0 ⎢ ⎢0  BI = ⎢ ⎢0 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎣

0 0

0 0

..

..

. ... ... ... ...

⎤

.

..

0 0 0 0

0 1 1 0

. 0

0 0 0 0

... ... ... ... .. .. . .

−1 −1 0

0 0

0

  u v

0 ⎢0

⎥ ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ ⎥, 0 ⎥ ⎥ 0 ⎥ .. ⎥ .⎥ ⎥ 0 ⎦

For more compact notation we also introduce

w=

⎡

⎤

S Br =



r2x + r2y Tr 0

BSr 0



0

B Sr

0 0

⎥ ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ ⎥. 0 ⎥ ⎥ 0 ⎥ .. ⎥ .⎥ ⎥ 0 ⎦ 0

 (50)

,

0 r2x + r2y Tr

⎞ ⎠,

= P



P 0

0

P

 ,

H=



Hqr

0

0

Hqr

 . (51)

The semi-discrete approximations of the systems (1), (2) with weak enforcements of the interface conditions (7) reads

d2



dt2



ρu T  w − τ N =P H−1 B I Tr w − γ N H−1 S TBr B I w − γN H−1 DCq B I w − τ0 H−1 B I w. ρv

(52)

In the right hand side of (52) the first term approximates the right hand sides of (1), (2), the second term enforces the continuity of stresses and the next three terms consistently enforce the continuity of displacements. The real numbers τ0 , τN , γN are penalty parameters determined by requiring stability. 4. Stability We will analyze the stability of the numerical interface treatment (52) for the 2nd, 4th and 6th order accurate SBP operators. The aim is to derive sharp estimates of the penalty terms τ0 , τ N , γ N such that we have a discrete energy estimate analogous to (32). To begin with, we set τ N = 12 and γ N = − 12 and use the SBP properties to rewrite (52) as

46

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

d2



dt2







1 −1 T 1 −1 T −D ρu 0 = H−1 w+ H  B I S Br w + H S Br BI w   ρv −D 0 2 2 1 −1 T 1 −1 T + H  B I DCq w + H DCq B I w − τ0 H−1 B I w, 2

where

(53)

2

    D = I 2 ⊗ D qT ⊗ I r ( I 2 ⊗ H q ⊗ H r ) A( I 2 ⊗ D q ⊗ I r ) + I 2 ⊗ I q ⊗ D rT ( I 2 ⊗ H q ⊗ H r )B( I 2 ⊗ I q ⊗ D r )     + I 2 ⊗ D qT ⊗ I r ( I 2 ⊗ H r ⊗ H q ) C( I 2 ⊗ I q ⊗ D r ) + I 2 ⊗ I q ⊗ D rT ( I 2 ⊗ H r ⊗ H q ) CT ( I 2 ⊗ D q ⊗ I r )



+ ( I 2 ⊗ I q ⊗ H r )Rq(A) + ( I 2 ⊗ H q ⊗ I r )Rr(B) ,     D  = I 2 ⊗ D qT ⊗ I r ( I 2 ⊗ H q ⊗ H r ) A ( I 2 ⊗ D q ⊗ I r ) + I 2 ⊗ I q ⊗ D rT ( I 2 ⊗ H q ⊗ H r ) B ( I 2 ⊗ I q ⊗ D r )     + I 2 ⊗ D qT ⊗ I r ( I 2 ⊗ H r ⊗ H q ) C ( I 2 ⊗ I q ⊗ D r ) + I 2 ⊗ I q ⊗ D rT ( I 2 ⊗ H r ⊗ H q ) C T ( I 2 ⊗ D q ⊗ I r ) 



+ ( I 2 ⊗ I q ⊗ H r )Rq( A ) + ( I 2 ⊗ H q ⊗ I r )Rr(B ) . By introducing

= D



D

0

0

D

 

+

1 2

(54)

(55)



T  2τ0 BI −  B TI S Br − S TBr BI −  B TI DCq − DCq BI ,

(56)

then the semi-discrete problem (53) becomes

d2 dt2





ρu w. = − H−1 D ρv

(57)

Note that the operator  B TI is defined by

⎧ if j = { N y , N y + 1},    ⎨0 u T  if j = N y , BI = ui , N y + vi ,1 ⎩ v ij −(ui , N y + vi ,1 ) if j = N y + 1.

(58)

The following theorem is central to the present study.

is symmetric positive semi-definite, then Theorem 1. Consider the system of ODEs (57). If D

     du 2  dv 2     w  + ρ  E w (t ) =  ρ + wT D dt dt  Hqr

(59)

Hqr

is a semi-discrete energy satisfying

E w (t ) = E w (0).

(60)

Proof. Multiply (57) with (dw/dt) T H from the left and add the transpose of the products we have

d dt

E w (t ) = 0.

(61)

Integrating in time yields the energy equation

E w (t ) = E w (0). The energy is conserved.

(62)

2

is symmetric. The difficulty lies in ensuring that the spatial operator D is positive semidefNote that by construction D , inite. We will show that the right choice of the penalty τ0 results to a symmetric positive semidefinite spatial operator D and thus yields an energy estimate. To begin with, we introduce the symmetric matrix

⎛

h Nr Nr B (qi , r Nr ) (r )

02

⎜ (r )  ⎜ 02 h11 B (qi , r1 ) ⎜ ⎜ (r ) L i = ⎜ h N N C (qi , r Nr ) 02 ⎜ r r ⎜ (r )  02 h11 C (qi , r1 ) ⎝ 1 1  B (q j , r Nr ) B (qi , r Nr ) 2 2

h Nr Nr C T (qi , r Nr )

02

02

h11 C  T (qi , r1 )

h Nr Nr A (qi , r Nr )

02

02

h11 A  (qi , r1 )

(r )

(r )

(r )

1 T C (qi , r Nr ) 2

(r )

1  T C (qi , r1 ) 2

⎞

1 B (qi , r Nr ) 2 ⎟ 1  B (qi , r1 ) ⎟ 2 ⎟ ⎟ 1 ⎟, C ( q , r ) Nr i 2 ⎟ ⎟ 1  C ( q , r ) ⎠ 1 i 2

τ0 I2

(63)

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Table 1 The constant

47

γ of Lemma 2.

Order

γ

2 4 6

0.5 0.3172 0.2950

where A, B, C, A , B, C  are the transformed coefficients matrices defined in (19) and (20). We will use the following lemma: Lemma 2. Consider the symmetric 10 × 10 matrix L i , defined in (63) above. There exists a constant

γ > 0, given in Table 1, such (b) (b ) (qi , r1 ))/(4γ h), (θ p (qi , r Nr ) + θ p (qi , r1 ))/(4γ h)), where h > 0 is the spatial step and the (b) (b ) eigenvalues θ , θ (with k = p , s) are defined in (23), (24), then the matrix L i is symmetric positive semi-definite. (b )

(b)

that if τ0 ≥ maxi ((θs (qi , r Nr ) + θs k

k

(r )

(r )

Proof. Consider h11 = h N r N r = γ h, where γ is given in Table 1 for SBP operators of order of accuracy 2, 4, 6. Note that by construction L i is symmetric. Consider now

U T L i U,

U = (U1 , U2 , U3 , U4 , U5 ) T ,

(64)

where U j = (U j1 , U j2 ) T , and j = 1, 2, . . . , 5. Introduce the nonsingular matrix

⎛

I2

⎜ 02 ⎜ ⎜ V i = ⎜ 02 ⎜ ⎝ 02 −1

02

02

I2

02

02

02

I2

02

02

02

I2

⎟ ⎟ 02 ⎠

I 2γ h 2

02

02

I2

−1

I 2γ h 2

02

⎞

02

02 ⎟ ⎟ 02 ⎟

(65)

and define

U = ViW

W = V i−1 U,

⇐⇒

we have

U T L i U = W T V iT L i V i W = W T L i W, and

⎛

γ h B

γ h CT

02

γ h B

02

γ h A

02

02

02

02

⎞

02

γ h C T

γ h A

(66)

i

02

⎜ 02 ⎜ ⎜ C 02 L i = ⎜ γ h ⎜ C ⎝ 02 γ h

02

Li = V T Li V i ,

⎟ ⎟ ⎟ ⎟. ⎟ ⎠

02

02

02 02

τ0 I2 −

B + B

(67)

4γ h

Thus, we have

UT L i U = WT Li W = γ h

 = γh

W3

T

W1

P



W3

T

W1



W3 W1

   T     B + B W3 W4 W4 +γh + W5T τ0 I2 − P P W5 W1



 +γh

W4 W2

W2

T

P



W4 W2

W2



⎛

+ ( S B W5 )T ⎝

4γ h

τ0 −

(b)

0

4γ h

0

⎞

(b )

θ p +θ p

τ0 −

(b)

(b )

θs +θs

⎠ ( S B W5 ).

4γ h

(b )

Here, θk , θk (with k = p , s) are the eigenvalues of B, B  and S B is the corresponding eigenvector matrix. Note that the (b) (b ) (b) potential energy matrices P, P  are symmetric positive semi-definite. Therefore, if τ0 ≥ max((θ p + θ p )/(4γ h), (θs +  (b ) (b)

)/(4γ h)), then the matrix L i is symmetric and positive semi-definite. Note that on rectangular Cartesian grids, the metric coefficients and the Jacobian are qq = 1, rq = 0, qr = 0, rr = 1 and J = 1. Thus, we have

θs

θ p(b) = c 22 ,



 θ p(b ) = c 22 ,

θs(b) = c 33 ,



 θs(b ) = c 33 .

2

48

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Introduce now the scalar products Nq Nr  

T

v ( H q ⊗ H r0 )u =

(q) (r )

h ii h j j v i j u i j ,

(68)

i =1 j =2

v T ( H q ⊗ H r N )u =

N q N r −1  

(q) (r )

h ii h j j v i j u i j ,

(69)

i =1 j =1

(q)

where we have excluded all points on the interface. Note that H q , H r are diagonal norms with H qii = h ii > 0, H r j j = (r )

h j j > 0. We define the quantity



wT Dw



interior

= uqT ( I 2 ⊗ H q ⊗ H r N ) Auq + uqT ( I 2 ⊗ H q ⊗ H r N ) Cur + urT ( I 2 ⊗ H q ⊗ H r N ) C T uq



( A )

( B )

(A) (B) + urT ( I 2 ⊗ H q ⊗ H r N ) Bur + u T ( I 2 ⊗ I q ⊗ H r )Rq u + u T ( I 2 ⊗ H q ⊗ I r )Rr u + vqT ( I 2 ⊗ H q ⊗ H r0 ) A vq + vqT ( I 2 ⊗ H q ⊗ H r0 ) C vr + vrT ( I 2 ⊗ H q ⊗ H r0 ) C T vq

+ vrT ( I 2 ⊗ H q ⊗ H r0 ) B vr + v T ( I 2 ⊗ I q ⊗ H r )Rq v + v T ( I 2 ⊗ H q ⊗ I r )Rr

v,

(70)

where we have excluded some of the interface terms in w T Dw. Dw)interior defined in (70) is nonnegative. Lemma 3. The quantity (w T Proof. Consider the expression (70) and simplify further, we have



T

w Dw

 interior

=

N q N r −1    uri j T i =1 j =1

uqi j

    Nq Nr    vri j T uri j vri j (q) (r ) (q) (r )  P (qi , r j ) h ii h j j + P (qi , r j ) h ii h j j uqi j

i =1 j =2



vqi j

vqi j





+ uT ( I 2 ⊗ I q ⊗ H r )Rq(A) u + uT ( I 2 ⊗ H q ⊗ I r )Rr(B) u + vT ( I 2 ⊗ I q ⊗ H r )Rq(A ) v ( B )

+ vT ( I 2 ⊗ H q ⊗ I r )Rr

v.

(71) ( A)

( B)

The matrices ( I 2 ⊗ H q ⊗ I r ) and ( I 2 ⊗ I q ⊗ H r ) are diagonal, and by (44), (45) the products ( I 2 ⊗ I q ⊗ H r )Rq , ( I 2 ⊗ H q ⊗ I r )Rr , ( A )

( B )

( A )

( B )

( I 2 ⊗ I q ⊗ H r )Rq , ( I 2 ⊗ H q ⊗ I r )Rr

( A)

( B)

commute. Thus, the corresponding matrices ( I 2 ⊗ I q ⊗ H r )Rq , ( I 2 ⊗ H q ⊗ I r )Rr ,

are symmetric positive semi-definite. Since the potential energy matrices P, P  are T also symmetric positive semi-definite, then it follows that (w Dw)interior ≥ 0. 2

( I 2 ⊗ I q ⊗ H r )Rq , ( I 2 ⊗ H q ⊗ I r )Rr

For convenience, we will omit the spatial arguments. The following theorem can be proven.

defined in (56). Let S r = D r , there exists a constant Theorem 2. Consider the matrix D

γ > 0, given in Table 1, such that if   is symmetric positive semi-definite, with τ0 ≥ max((θ p(b) + θ p(b ) )/(4γ h), (θs(b) + θs(b ) )/(4γ h)), where h > 0 is the spatial step, D wT Dw ≥ 0 for all w ∈ Rn .

defined in (56) is symmetric. To complete the proof of the lemma, we will Proof. By construction the spatial operator D is positive semi-definite. To begin, replace S r with D r in (56), and introduce prove that D

 

uq = ( I 2 ⊗ D q ⊗ I r )u,

ur = ( I 2 ⊗ I q ⊗ D r )u,

w=

u v

.

Consider

wT Dw = uqT ( I 2 ⊗ H q ⊗ H r ) Auq + uqT ( I 2 ⊗ H q ⊗ H r ) Cur + urT ( I 2 ⊗ H q ⊗ H r ) C T uq ( A) ( B) + urT ( I 2 ⊗ H q ⊗ H r ) Bur + u T ( I 2 ⊗ I q ⊗ H r )Rq u + u T ( I 2 ⊗ H q ⊗ I r )Rr u + vqT ( I 2 ⊗ H q ⊗ H r ) A vq + vqT ( I 2 ⊗ H q ⊗ H r ) C vr + vrT ( I 2 ⊗ H q ⊗ H r ) C T vq 



(A ) (B ) + vrT ( I 2 ⊗ H q ⊗ H r ) B vr + v T ( I 2 ⊗ I q ⊗ H r )Rq v + v T ( I 2 ⊗ H q ⊗ I r )Rr v Nq

+

 (q) 1  2τ0 (vi1 − uiN r ) T (vi1 − uiN r ) + (vi1 − uiN r ) T BuriN r + (vi1 − uiN r ) T B vri1 h ii 2 i

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Nq

+

1  2

 (q) (vi1 − uiNr )T C T uqiN r + (vi1 − uiN r ) T C T vqi1 h ii

i Nq

+

1  2



(q)



(q)

T T  uriN B(vi1 − uiN r ) + vri1 B (vi1 − uiN r ) h ii r

i Nq

+

49

1  2

T T  uqiN C(vi1 − uiN r ) + vqi1 C (vi1 − uiN r ) h ii . r

i

Introducing the scalar products (68) and (69), where we have excluded all points on the interface, we have

⎛

uriN r

Nq ⎜  ⎜ ⎜ w Dw = w Dw interior + ⎜ i ⎝



T

T

vri1



uqiN r vqi1

⎞T

⎟ ⎜ ⎟ ⎜ ⎟ Li ⎜ ⎟ ⎜ ⎠ ⎝

(vi1 − uiNr ) (b )

(b)

(b)

⎛

uriN r vri1 uqiN r vqi1

⎞ ⎟ ⎟ (q) ⎟h . ⎟ ii ⎠

(72)

(vi1 − uiNr )

(b )

If τ0 ≥ max((θ p + θ p )/(4γ h), (θs + θs )/(4γ h)), where h > 0 is the spatial step, then the matrix L i is symmetric positive semi-definite. Therefore, the interior terms and the interface terms in (72) are nonnegative, and w T Dw ≥ 0. The proof of the theorem is complete. 2 Since we have assumed S r = D r , the result above holds for fully compatible SBP operators. To be able to extend the analysis to more general classes of compatible SBP operators with S r = D r we introduce (a)

(a)

R sr = R r

  T  − ( E Rr − E Lr )( S r − D r ) H r a ( E Rr − E Lr )( S r − D r ) ,

(73)

(a)

(a)

where R r denotes the corresponding remainder operator of a minimal width stencil SBP operator D rr and E Rr = diag(0, 0, . . . , 1), E Lr = diag(1, 0, . . . , 0). For the special and important case of constant coefficients (a ≡ 1) the operators (a) (a) R r , R sr have simple expressions. We present below the corresponding operators for the second order accurate SBP operator,

⎡

1 −2 ⎢ −2 5 ⎢ ⎢ 1 −4

(a)

Rr

=

1 ⎢

⎢

..

4h ⎢ ⎢

.

⎢ ⎣ ⎡

1 −4 1 6 −4

..

.

1

1 2

−1

1 2

⎤

⎥ ⎥ ⎥ ⎥ .. .. .. ⎥, . . . ⎥ ⎥ −4 6 −4 1 ⎥ ⎦ 1 −4 5 −2 1 −2 1 1

0

⎤

⎢ −1 3 −3 1 ⎥ 0 ⎢ ⎥ ⎢ 1 −3 11 −4 1 ⎥ 0 ⎢ 2 ⎥ 2 ⎢ ⎥ ⎢ 0 ⎥ 0 1 −4 6 −4 1 ⎢ ⎥ 1 ⎢ ⎥ (a) . . . . . . . .. .. .. .. .. .. .. R sr = ⎢ ⎥. ⎥ 4h ⎢ ⎢ 0 1 −4 6 −4 1 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ −3 12 ⎥ 0 1 −4 11 2 ⎢ ⎥ ⎣ 0 1 −3 3 −1 ⎦ 1 −1 12 0 2 The following theorem was proven in [22]. n

Theorem 3. (See Theorem 3.3 in [22].) Let A be an n × n penta-diagonal symmetric matrix. Assume j =1 A i j = 0 and A ii > 0, A i ,i +1 < 0, A i ,i +2 > 0. If − A 1,2 ≥ 2 A 1,3 , − A n,n−1 ≥ 2 A n,n−2 and − A i ,i +1 ≥ 2 A i −1,i +1 + 2 A i ,i +2 , i = 2, . . . , n − 2, then A is symmetric positive semi-definite. (a)

Theorem 3 above was used in [22] to prove the positive semi-definiteness of the remainder operator R r , for the second (a) (a) order accurate case. Note that both R r and R sr , above, satisfy all the conditions in Theorem 3. Therefore, the operator

50

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

(a)

R sr is also symmetric positive semi-definite. The corresponding results for the fourth and sixth order accurate schemes in (a) [22] can also be extended to R sr . A direct consequence of Theorem 3 above and Proposition 3.5 in [22] yield the following lemma (a)

Lemma 4. Consider the minimal width stencil SBP operator D rr for the second derivative ∂/∂ r (a∂/∂ r ), defined in (34), with S r = D r (a) being the discrete approximation of the first derivative ∂/∂ r on the boundaries. Let R r denote the corresponding remainder operator and E Rr = diag(0, 0, . . . , 1), E Lr = diag(1, 0, . . . , 0) are boundary matrices for the right and left boundaries respectively. The operator (a) R sr defined by (73) is symmetric and positive semi-definite. For convenience, we will omit the spatial arguments (xi , y N y ). We will also use the symmetric matrix

⎛

K=

2τ02 I 2

1⎜

B

⎝

2

B

B

2h N r N r B

0

(r )

B

⎞

(r ) 

0

⎟ ⎠,

(74)

2h11 B

and the corresponding lemma Lemma 5. Consider the 6 × 6 matrix K defined in (74) above. There exists a constant (b)

(b )

(b)

(b )

max(θ p /(2γ h), θs /(2γ h), θ p /(2γ h), θs definite.

γ > 0, given in Table 1, such that if τ02 ≥ /(2γ h)), where h > 0 is the spatial step, then the matrix K is symmetric positive

Proof. Consider

T



p q r

p q r

K

=

 T

1

 

p q

2

K1

p q

+

where K 1 , K 2 are 4 × 4 matrices defined by



K1 =

B



B

τ02 I 2

(r )

2h N r N r

K2 =

, B

1

 T p r

2



τ02 I 2 B

  p r

K2

,



B

(75)

.

2h11 B (r )

We can simplify further having

 T

 

p q

and

p q

K1

 T

q1

 

p r

p r

K2

  T      T    p1 p1 p2 p2 = SB + SB , K 11 S B K 12 S B q1

K 11 =

τ02 (b)

K 22 =

θs

τ02 (b )

θp

q2

  T      T    p1 p1 p2 p2 = SB + SB . K 21 S B K 22 S B r1

r1

Here, S B are the eigenvector matrices and



q2

(r )





(b)

θs

,

(b)

2h N r N r θs 

θ p(b ) (r ) (b )

K 12 =

τ02 (b)

θp



r2





(b)

θp (r )

r2

(b)

2h N r N r θ p

,

K 21 =

τ02

(b )

θs



λ

K 12

=

τ02 + 2θs(b) h(Nrr)Nr ± (τ02 + 2θs(b) h(Nrr)Nr )2 + 4((θs(b) )2 − 2τ02 θs(b) h(Nrr)Nr ) τ02 + 2θ p(b) h(Nrr)Nr ± (τ02 + 2θ p(b) h(Nrr)Nr )2 + 4((θ p(b) )2 − 2τ02 θ p(b) h(Nrr)Nr ) 

λ K 21 = λ

=

2 





τ02 + 2θs(b ) h(11r ) ± (τ02 + 2θs(b ) h(11r ) )2 + 4((θs(b ) )2 − 2τ02 θs(b ) h(11r ) ) 

K 22

2





(r ) (b )

2h11 θs

 , (76)

,

2h11 θ p

are 2 × 2 matrices. The eigenvalues of the matrices K 11 , K 12 , K 21 , K 22 are

λ K 11 =

(b )

θs

2 





τ02 + 2θ p(b ) h(11r ) ± (τ02 + 2θ p(b ) h(11r ) )2 + 4((θ p(b ) )2 − 2τ02 θ p(b ) h(11r ) ) 2

, .

, ,

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

(r )

(r )

(b)

(b )

(b)

(b )

Therefore, if h11 = h N r N r = γ h and τ02 ≥ max(θ p /(2γ h), θs /(2γ h), θ p /(2γ h), θs K 21 , K 22 are symmetric and positive definite. Thus, implying

T p q r

51

/(2γ h)) then the matrices K 11 , K 12 ,



K

p q r

≥ 0,

for any p, q, r ∈ R2 .

2

The following and more general result can be proven

defined in (56). For S r = D r , if Theorem 4. Consider the matrix D























τ0 ≥ max θ p(b) + θ p(b ) /(4γ h), θs(b) + θs(b ) /(4γ h) + max θ p(b) /(2γ h), θs(b) /(2γ h), θ p(b ) /(2γ h), θs(b ) /(2γ h) , is symmetric positive semi-definite, with w T where h > 0 is the spatial step and γ given in Table 1, then the spatial operator D Dw ≥ 0 for all w ∈ Rn . defined in (56) is symmetric. It remains to prove that D is positive definite. Proof. We know that the spatial operator D Introduce us = ( I 2 ⊗ I q ⊗ S r )u,

ur = ( I 2 ⊗ I q ⊗ D r )u,





usr = I 2 ⊗ I q ⊗ ( S r − D r ) u,

with us = ur + usr . Consider

wT Dw = uqT ( I 2 ⊗ H q ⊗ H r )Auq + uqT ( I 2 ⊗ H q ⊗ H r )Cur + urT ( I 2 ⊗ H q ⊗ H r )C T uq

+ urT ( I 2 ⊗ H q ⊗ H r )Bur + uT ( I 2 ⊗ I q ⊗ H r )Rq( A ) u + uT ( I 2 ⊗ H q ⊗ I r )Rr( B ) u + vqT ( I 2 ⊗ H q ⊗ H r )A vq + vqT ( I 2 ⊗ H q ⊗ H r )C vr + vrT ( I 2 ⊗ H q ⊗ H r )C T vq 



+ vrT ( I 2 ⊗ H q ⊗ H r )B vr + vT ( I 2 ⊗ I q ⊗ H r )Rq( A ) v + vT ( I 2 ⊗ H q ⊗ I r )Rr( B ) v Nq

+

1  2



2τ0 (vi1 − uiN r ) T (vi1 − uiN r ) + (vi1 − uiN r ) T BusiN r + (vi1 − uiN r ) T B  vsi1 h ii

i Nq

+

1  2

 (x) (vi1 − uiNr )T C T uxiNr + (vi1 − uiNr )T C  T vxi1 hii

i Nq

+

1  2

(q)



(q)

i

1  2



T T usiN B (vi1 − uiN r ) + vsi1 B  (vi1 − uiN r ) h ii r

Nq

+

(x)

T T uqiN C (vi1 − uiN r ) + vqi1 C  (vi1 − uiN r ) h ii . r

i

Introducing the scalar products

v T ( H q ⊗ H r0 )u =

Nq Nr  

(x) ( y )

h ii h j j v i j u i j ,

v T ( H q ⊗ H r N )u =

i =1 j =2

N q N r −1  

(x) ( y )

h ii h j j v i j u i j ,

i =1 j =1

where we have excluded all points on the interface, we have

wT Dw = uqT ( I 2 ⊗ H q ⊗ H r N )Auq + uqT ( I 2 ⊗ H q ⊗ H r N )Cur + urT ( I 2 ⊗ H q ⊗ H r N )C T uq

+ urT ( I 2 ⊗ H q ⊗ H r N )Bur + uT ( I 2 ⊗ I q ⊗ H r )Rq( A ) u + uT ( I 2 ⊗ H q ⊗ I r )R(srB ) u + vqT ( I 2 ⊗ H q ⊗ H r0 )A vq + vqT ( I 2 ⊗ H q ⊗ H r0 )C vr + vrT ( I 2 ⊗ H q ⊗ H r0 )C T vq 



+ vrT ( I 2 ⊗ H q ⊗ H r0 )B vr + vT ( I 2 ⊗ I q ⊗ H r )Rq( A ) v + vT ( I 2 ⊗ H q ⊗ I r )R(srB ) v ⎞T ⎛ ⎞ ⎛ uriN r uriN r ⎛ ⎞ ⎛ ⎞ ⎟ ⎟ ⎜ Nq ⎜ Nq vri1 vri1 (vi1 − uiNr ) T (vi1 − uiNr )  ⎟ ⎟ (q)  ⎜ ⎜ ⎟ L i ⎜ uqiN ⎟h + ⎜ uqiN ⎝ usriNr ⎠ K ⎝ usriNr ⎠ h(q) . + r r ii ⎟ ⎟ ii ⎜ ⎜ ⎠ ⎠ ⎝ ⎝ i i v v v v qi1

(vi1 − uiNr )

qi1

(vi1 − uiNr )

sri1

sri1

(77)

52

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Note we have used the identity S r = D r + ( S r − D r ) ⇒ us = ur + usr and (B)

Rr

= R(srB ) +





T

I 2 ⊗ I q ⊗ ( E R y − E L y ) (Sr − Dr )

   ( I 2 ⊗ I q ⊗ H r )B I 2 ⊗ I q ⊗ ( E R y − E L y ) (Sr − Dr ) .

Thus, if























τ0 ≥ max θ p(b) + θ p(b ) /(4γ h), θs(b) + θs(b ) /(4γ h) + max θ p(b) /(2γ h), θs(b) /(2γ h), θ p(b ) /(2γ h), θs(b ) /(2γ h) , where h > 0 is the spatial step, then the matrices L and K are symmetric positive semi-definite. By Lemma 4 we have that (B) (B) Rsr = (Rsr ) T ≥ 0. Therefore, the interior terms and the interface terms in (77) are nonnegative, yielding w T Dw ≥ 0. On rectangular Cartesian grids, the metric coefficients and the Jacobian are qq = 1, rq = 0, qr = 0, rr = 1 and J = 1. Thus, we also have 

 θ p(b ) = c 22 ,

θ p(b) = c 22 ,



 θs(b ) = c 33 .

θs(b) = c 33 ,

2

Note that when S r = D r we need a larger penalty strength, τ0 , thus increasing the stiffness of the system of ODEs (53). Although, the fully-compatible SBP operator (with S = D r ) reduces the order of accuracy of the numerical approximation by one, at the boundary points, the stability analysis of the corresponding scheme is less arduous as compared to the compatible SBP operator (with S = D r ). Interestingly, numerical experiments performed in the next section show that the reduction in the accuracy of the boundary stencil does not affect the global accuracy of the scheme. When the compatible SBP operator is used, the order of accuracy on the boundary stencil is preserved. However, the analysis above shows that a slightly larger penalty strength is needed for stability. We have also performed numerical experiments using the larger penalty parameters derived above for the compatible SBP operators. Numerical experiments show that the slight increase in the penalty strength does not restrict the time step of an explicit scheme further. 5. Numerical experiments In this section, we present numerical experiments to verify the theory developed in this paper. We use schemes that have been constructed using 4th and 6th order diagonal fully compatible and compatible SBP operators. The material parameters are chosen to represent isotropic elastic media i.e., if λ and μ are the first and second Lamé parameters, respectively, then

c 11 = c 22 = λ + 2μ,

c 12 = λ,

c 33 = μ.

Note that in curvilinear coordinates the corresponding coefficient matrices have a more general form. To integrate numerically in time we use the 4th order scheme described in [8,10] which was designed for second order systems of ODEs,

U tt = MU resulting from semi-discretization of second order wave equations using the SBP–SAT scheme. It can be proved [10] that if the time step t satisfies

t ≤

2.54

(M1 M∞ )1/4

(78)

,

then the schemes are stable. In the experiments, unless otherwise stated, equality was chosen in (78). We present five numerical experiments. The first two serve as verifying accuracy and stability. The third experiment displays the capability of handling simulations that includes the Stoneley interface wave, the fourth studies artificial reflections at domain interfaces and the last experiment demonstrates the geometric capabilities of the proposed method. 5.1. Mode conversion at a line interface In the first experiment we consider a compressional plane wave of unit amplitude propagating with angle θ and temporal frequency 2ωπ in the negative y-direction. The displacement field is given by



(in)

UP

=

 ξ e i (γ ξ x−γ η y −ωt ) , −η

  π ω , ξ = sin(θ), η = cos(θ), γ = √ θ ∈ 0, , x ∈ (−∞, ∞), y ∈ [0, ∞). 2 λ + 2μ

We assume an interface between two different materials at y = 0. The Lame parameters and densities are λ, μ, ρ and λ , μ , ρ  in the half-planes y > 0 and y < 0, respectively. When this wave encounters the interface between the two different materials it will be split into reflected, compressional and shear waves,

U =U

(in)

(refl)

+ UP

η2 = cos(θ2 ),

(refl)

+ US

,

sin(θ2 ) =

(refl)

UP

γ ξ, γ2

= A refl

  ξ

η

e

i (γ ξ x+γ η y −ωt )



,

(refl)

US

ω γ2 = √ , x ∈ (−∞, ∞), y ∈ [0, ∞), μ

= B refl



γ2 η2 i(γ ξ x+γ2 η2 y−ωt ) e , γξ

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

53

Fig. 1. Max error as a function of time for a plane wave impinging on an interface where discontinuities in the Lamé parameters occur. Results from methods using fourth and sixth order SBP operators are shown from left to right. The number of grid points increases from top to bottom.

and refracted compressional and shear waves,

  γ1 ξ1 i(γ1 ξ1 x−γ3 η3 y−ωt ) (refr) (refr) + US , UP = A refr e , γ 3 η3   −γ4 η4 i (γ1 ξ1 x−γ4 η4 y −ωt ) (refr) US = B refr e , γ 1 ξ1 U = UP

(refr)

η3 = cos(θ3 ),

sin(θ3 ) =

x ∈ (−∞, ∞),

γ1 ξ1 , γ3

γ3 = 

ω λ

+ 2μ

,

η4 = cos(θ4 ),

sin(θ4 ) =

γ1 ξ1 , γ4

ω γ4 =  , μ

y ∈ (−∞, 0].

The constants A refl , B refl , A refr , B refr are obtained by inserting U and U  into the conditions (7) and solving the resulting linear system. For a more detailed discussion see [11, pp. 377–380]. We choose θ = π4 , ω = 2π , λ = μ = 1, λ = 1/2, μ = 2 and ρ = ρ  = 1. We use this exact solution to verify stability and accuracy of the numerical schemes described in this work. The computational domain is taken to be (x, y ) ∈ [− 2γπξ , 2γπξ ] × [− γ2πη , γ2πη ]. Initial data for the numerical scheme is 2 2

3 3

taken as the real part of the analytic solution at time t = 0 and exact data is imposed at the outer boundaries. The solution is computed for 100 periods until t = 100 and the discrete maximum error is measured. We use 2N × N grid points. Figs. 1(a)–1(b) display the maximum error in the numerical solution as a function of time for t ≤ 100 obtained with schemes using compatible SBP operators. Figs. 1(c)–1(d) display the corresponding measurements obtained with fully compatible SBP operators. These figures illustrate the stability and accuracy of the schemes. Numerical values of the rate of convergence measured at time t = 100 are shown in Table 2. Note that the computations for this example suggest that the 6th order fully compatible SBP operators yield schemes that have slightly better accuracy and convergence properties than the corresponding 6th order compatible versions even though the formal order of accuracy is one order higher at exactly one point for the compatible SBP operators. 5.2. A curvilinear interface In this section, we present numerical experiments for a two-block elastic medium separated by a curved interface. We will demonstrate the applicability of our method to curved interfaces and numerically evaluate accuracy and stability.

54

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Table 2 Rate of convergence of the solution in max norm. p 4 and p 6 give the measured convergence using 4th and 6th order operators, respectively. Compatible SBP operators

Fully compatible SBP operators

N

p4

p6

p4

p6

21 41 81 161

– 3.98 4.01 4.00

– 6.53 5.57 5.52

– 4.07 4.00 3.99

– 7.06 5.98 6.37

Fig. 2. A two-block curvilinear grid and a plot of the point wise error on the grid. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

To begin with, we consider two elastic blocks separated by an interface at

y = 1 + 0.2 sin(2π x).

(79)

We split the block into two, a lower block and an upper block, and generate a conforming curvilinear grid on each block. A plot of the grid using 21 grid-points on each block is shown in Fig. 2(a). Note that the curved interface is painted red. We map the lower and upper curvilinear blocks separately to a computational unit square. We have imposed traction boundary conditions on the outer boundaries. In order to quantify numerical errors, we have chosen forcing functions, the initial conditions and boundary data to satisfy the exact solution

u = cos(2π t ) sin(2π x) sin(2π y ),

v = − cos(2π t ) sin(2π x) sin(2π y ).

(80)

We discretize the equations using compatible and fully compatible SBP operators. The compatible SBP operators for variable coefficients are derived in [17]. To the best of our knowledge, the fully compatible SBP operators for variable coefficients have not been presented in the literature. However, as shown in Lemma 1 on p. 43 of this paper, these operators can be derived by modifying the boundary derivative operators corresponding to the compatible SBP operators derived in [17]. We have used SBP operators of interior order of accuracy 4 and 6, respectively. The interface conditions and all boundary conditions are imposed weakly using SAT. The solution is advanced in time using a conservative fourth order accurate time-stepping scheme [8]. We consider an isotropic elastic medium and use λ = 2, μ = 1, ρ = 1 in both blocks. In Fig. 2(b), we have plotted a snapshot (at t = 0.1) of the point wise error on the grid using a sixth order accurate SBP operator with the spatial step h = 0.0125. Note that the maximum errors on the grid occur where the one-sided stencils corresponding to the first derivative operators and their transpose interact in a non-trivial way. We have also performed convergence tests. In Figs. 3 and 4, we display the numerical errors for the compatible and fully compatible operators, respectively, at different resolutions. Note that numerical errors obtained using the compatible and fully compatible operators are analogous and converge to zero at the theoretical convergence rate predicted in [28] (rate = min(r + 2, 2r ), with r = 2, 3 for the fourth order and sixth order accurate schemes, respectively). 5.3. The Stoneley interface wave A Stoneley interface wave clings to the interface between two elastic media and decays exponentially, into the media, away from the interface. Define the p-wave speed and s-wave speed

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

55

Fig. 3. Time series of the maximum error showing the convergence of global numerical errors arising from the numerical treatment of a curvilinear interface using compatible SBP operators.

Fig. 4. Time series of the maximum error showing the convergence of global numerical errors arising from the numerical treatment of a curvilinear interface using fully compatible SBP operators.

!

α=

λ + 2μ

ρ

" ,

β=

μ . ρ

To simplify, consider two elastic half-planes with an interface at y = 0 and assume that the Stoneley wave is 2π -periodic in the x-direction. The component of the Stoneley wave in the positive half-plane y ≥ 0 can then be written

U = Ae

− y 1−c 2S /α 2





cos(x − c S t )

1 − c 2S /α 2 sin(x − c S t )



+ Be

− y 1−c 2S /β 2

 − 1 − c 2S /β 2 cos(x − c S t ) − sin(x − c S t )

and the component of the Stoneley wave in the negative half-plane y ≤ 0 is written



U = Ce

y 1−c 2S /α  2



cos(x − c S t )

− 1 − c 2S /α  2 sin(x − c S t )





+ De

y 1−c 2S /β  2



1 − c 2S /β  2 cos(x − c S t )

− sin(x − c S t )

y ≥ 0,

,

(81)

 ,

y ≤ 0.

(82)

As in the experiment in Section 5.1 the constants A, B, C and D are determined through the solution to the linear system arising from inserting (81) and (82) into the interface conditions (7). Let c S denote the phase velocity of the Stoneley wave. Solutions of the form (81) and (82) exist if and only if c S satisfies the dispersion relation

56

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Fig. 5. The Stoneley surface wave as a function of (x, y ) at t = 0. The u 1 component is show to the left and the u 2 component to the right.

# # # 1 # # " # # c2 1 − αS2 # # det # # # −ρ c 2 − 2ρ β 2 # S # " # 2 # # −2ρ β 2 1 − c S2 α

" − 1−

β2

−1 " 2ρ β

2

"

c 2S

"

c 2S

1 − β2 c 2S

ρ  c2 S

c2

+ 2ρ β  2

−2ρ  β  2

"

# # # # # # # 1 # # " # = 0. 2 # cS #   2 2ρ β 1 − β 2 # # # 2 # c   2 S −ρ β (2 −  2 ) # − 1−

1 − α S2

ρ β (2 − β 2 ) 2

−1

c2

1 − α S2

c 2S

β 2

(83)

β

The existence of a root c S to Eq. (83) was first considered by Stoneley [27], where existence was shown for some special values of the choice of Lamé parameters and densities of the half-planes. In the following experiments a root will exist for all media parameters considered. The solutions (81) and (82) now represent waves traveling harmonically along the interface between the half-planes and decays exponentially into the domain. The initial experiment is aimed at evaluating the convergence properties of the current method when simulating Stoneley waves. In this experiment we will also study the difference in accuracy between schemes using 6th order compatible and fully compatible operators observed in the experiment in Section 5.1. To begin with, we take λ = μ = ρ = 1 and λ = 3, μ = 2, ρ  = 1.98. For this case a root

c S = 0.999240140585103 to (83) is found. Fig. 5 displays vertical and horizontal components of the Stoneley interface wave for this choice of parameters. The computational domain is taken to contain exactly one wavelength, 2π , in the x-direction. Periodic boundary conditions are imposed at x = 0 and x = 2π and exact Dirichlet data are imposed at y = 4π and y = −4π . Initial data for the computations is taken from (81) and (82) at time t = 0. We use a grid size of h = 2π /( N − 1) for N = 21, 41, 81, 161 and the max error in the approximate solutions are computed as a function of time for 5 temporal periods. The temporal period is given by

2π /c S ≈ 6.29. The computations are done for each N and each of the two schemes. Fig. 6 compares the errors obtained with each scheme. A dashed and full lines correspond to the scheme using fully compatible SBP operators and compatible SBP operators, respectively. It becomes evident that the scheme using fully compatible operators converges faster and has a smaller error. For N = 161 the error is more than one magnitude smaller. The Rayleigh wave is similar to the Stoneley interface wave, in fact it can be considered a special case of the Stoneley interface wave when one of the half-planes is a vacuum. It has been shown by [16] that in numerical simulations involving the Rayleigh wave an unexpectedly large number of grid points per wavelength is needed for accuracy. In [16] it is shown that the number of grid points per wavelength of the Rayleigh wave scales as (λ/μ)1/ p , where p is the order of accuracy of the method. For materials with λ  μ this scaling becomes large. The purpose of this experiment is to evaluate the performance of the current method when simulating Stoneley waves in materials with λ  μ and to further study the difference in accuracy between schemes using 6th order compatible and fully compatible. Here, we take λ = ρ = 1, λ = 3, λ = 2, ρ  = 1.9999, μ = 0.1, 0.01, 0.001, 0.0001 and μ = 2μ. For all considered values of parameters the exponential decay is of the same order of magnitude.

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

57

Fig. 6. Max errors as functions of time. A dashed line corresponds to the errors obtained with a scheme using fully compatible 6th order operators. A full line corresponds to the errors obtained with a scheme using compatible 6th order operators.

Fig. 7. Maximum errors as functions of time. Note that the grid size is the same for all computations. The magnitude of the errors are more than one magnitude smaller using fully compatible SBP operators (right) than those obtained using compatible SBP operators (left).

As above, the computational domain is taken to contain exactly one wavelength, 2π , in the x-direction. Periodic boundary conditions are also imposed at x = 0 and x = 2π and exact Dirichlet data are imposed at y = 4π and y = −4π . We let the wave be well resolved in the spatial coordinates by choosing a grid size of h = 2π /(80) this amounts to 81 grid points per wavelength. The initial data for the computations are taken from (81) and (82) at time t = 0 and we compute for 2 periods in time. Figs. 7(a)–7(b) display the maximum errors as a function of time. We see that as the quotients λ/μ and λ /μ become larger the magnitude of the error increases even though the waves are equally well resolved. That is, in simulations with λ  μ the Stoneley wave as the Rayleigh wave also requires an unexpectedly large resolution for an accurate solution. From the scaling (λ/μ)1/ p , the amplification of the magnitude of the error is worse for lower order schemes, see [7,16]. We conclude that high order accurate methods are essential for efficient simulations of Stoneley waves, for material with high Poisson ratio, when λ  μ. Further, by comparing the magnitude of the maximum errors between the two different schemes we see that the errors obtained using the fully compatible SBP operators are about one magnitude smaller than those obtained using the compatible SBP operators.

58

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Fig. 8. The horizontal component of the planar wave and the computational domain used in Experiment 4. The red lines illustrate how the domain is composed into four blocks and then patched together. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5.4. Artificial reflections It is important that artificial reflections arising from the numerical treatment of interface conditions at sub-domain interfaces do not contaminate the numerical solution. To measure such reflections we consider a planar pressure wave traveling in the plane. The (u , v ) components of the wave are expressed as 2

u=−

cos θ − (cos θ x+sin θ 2y−1−c p t ) 2σ e (cos θ x + sin θ y − 1 − c p t ), 2

σ

2

sin θ − (cos θ x+sin θ 2y−1−c p t ) 2σ v =− 2 e (cos θ x + sin θ y − 1 − c p t ).

σ

√

√

(84)

√

We take λ = μ = ρ = 1, θ = π /4 and σ = 1/300. Then c p = (λ + 2μ)/ρ = 3. With f 0 = 2π1σ ≈ 2.7566 we approximate the fundamental frequency of the wave as F = 2.5 f 0 ≈ 6.8916. The dominant wavelength of the wave is then Λ = c p / F ≈ 0.2513 and its period is 1/ F ≈ 0.1451. We let the computational domain be the square [0, 2] × [0, 2], the wave and the computational domain is displayed in Fig. 8. To measure reflections from the interface treatment described in this paper we cut the computational domain into four squares of unit length. We discretize the equations in the subdomain independently and then patch the subdomains into the global domain using penalties. The wave is propagated for four periods in time. At this time the wave has traveled across all four interfaces in the domain. We use both 4th order compatible and fully compatible SBP operators in this experiment. To give a quantitative statement of the size of the reflections we follow the classical theory concerning accuracy and grid points per wavelength when simulating hyperbolic wave propagation [15]. In particular, for a 10% and 1% relative max error using a 4th order method it is claimed that 7q1/4 and 13q1/4 grid points per wavelength, respectively, should be used. Here q is the number of periods the wave is being propagated. We therefore propose that in order for the interface treatment to be useful, reflections should not be greater than 10% or 1% of the magnitude of the true solution, corresponding to using 7q1/4 or 13q1/4 grid points per wavelength. We measure the relative max error as a function of time in the subdomain [0, 1] × [0, 1]. This part of the domain contains artificial reflections coming from interfaces to all other subdomains, including the corner points where the four sub-domains meet. Fig. 9 displays the errors, in both cases the reflections measured with the schemes using compatible SBP operators are less than the required error tolerance. The scheme using fully compatible SBP operators fails at fulfilling the 10% error tolerance but satisfies the 1% error tolerance. The failure in the first case is probably due to the accuracy reduction of the fully compatible operators. A time series of the error in the numerical solution obtained with the scheme constructed with fully compatible SBP operators and 7q1/4 grid-points per wavelength is displayed in Fig. 10. Here, the plot is zoomed in the region [1, 3] × [1, 3] to give a better view of the errors. 5.5. Scattering from a cavity with an irregular surface The propagation of a pressure wave in an elastic material can be described by the elastic wave equation. Numerical methods for scattering problems based on acoustic modeling are sometimes used [5]. However, as the pressure wave encounters a free surface parts of its energy is being transformed into shear waves, a phenomena not captured by the acoustic model. The boundaries of cavities in the domain of interest are in general non-planar. At cavity boundaries a traction free boundary condition applies. The purpose of this experiment is twofold. One: to illustrate the importance of using elastic

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

59

Fig. 9. Relative max errors obtained with the schemes using 4th order compatible and fully compatible SBP operators. The straight red line indicates the error tolerance. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Time series of the magnitude of the error showing the artificial reflections arising from the interface treatment.

modeling in simulations including scattering problems. Two: to display the geometric flexibility of the proposed method. This is achieved by considering a plane pressure wave impinging on an irregular shaped cavity, left of Fig. 11. The problem is solved with both an acoustic solver [31] and the elastic solver proposed in this paper. The solutions are presented by snapshots of the magnitude of the displacement field at three different times, Fig. 12. The computational domain is the box [−3, 3] × [−3, 9] with the cavity centered at the origin. The boundary ∂Ω of the cavity is described by





x = 1 + 0.2 cos(4π θ − π ) cos(π /2θ − π /4),





y = 1 + 0.2 cos(4π θ − π ) sin(π /2θ − π /4),

θ ∈ [0, 2π ).

The elastic material is defined by √ the Lamé parameters λ = 2, μ = 1 and time is scaled to give unit density. This choice gives a pressure wave speed c p = λ + 2μ = 2. Using the acoustic model the displacement field (u 1 , u 2 ) is computed from the solution to the scalar wave equation

φtt = c p ∇ 2 φ, ∂n φ = 0,

(x, y ) ∈ ∂Ω

u 1 = φx ,

(85)

u2 = φ y ,

(86)

where ∂n denotes the normal derivative. The incoming plane wave is initiated in the acoustic model by introducing a forcing at the boundary y = 9

  2 φ(x, 9) = f = sin 2π f c (9 + c p t ) e −2(c p t −5) .

(87)

Here, the central frequency f c is chosen to be f c = 2. In the elastic model the incoming plane wave is initiated by the boundary conditions

u 1 (x, 9) = 0

      2 u 2 (x, 9) = 2π f c cos 2π f c (9 + c p t ) − 4 sin 2π f c (9 + c p t ) (c p t − 5) e −2(c p t −5) ;

(88)

By (85)–(86) the solutions of both models are identical to machine precision sufficiently far away from the cavity. At the vertical boundaries the numerical solutions are enforced to coincide with the incoming plane wave, at the horizontal boundary

60

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

Fig. 11. The incoming plane wave and the computational domain. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

y = −3 the solutions are set to zero. To handle the geometry of the cavity the computational domain is decomposed into five blocks and then patched together with the SAT technique, Fig. 11. Each of the five blocks are resolved using 161 × 161 grid points and both the acoustic and the elastic solver use the same grid in the computations. In the computations fourth order discretizations are used. The time step is in both cases t = 5.5E−4. Snapshots of the magnitude of the computed displacement field are presented for both solutions at times t = 6.4, t = 7.0 and t = 7.6, Fig. 12. Both solutions correspond in location and frequency of the reflected pressure waves. In the vicinity of the cavity reflected shear waves of higher frequency are present in the solution obtained with the solver for the elastic wave equation described in this paper. This key feature in the reflected wave field is absent in the solution obtained with the acoustic solver. This observation is particularly visible in the snapshots of the solutions at time t = 7.6. 6. Conclusions High order accurate finite difference schemes for the two dimensional elastic wave equation in discontinuous and layered elastic media are presented. The key ingredient in the discretization was fully compatible or compatible finite-difference operators approximating first and second derivatives. Stability of the numerical scheme was proven for the discretization using energy methods. The convergence and accuracy properties, for the compatible and fully compatible operators, are observed to be similar in many cases. The numerical experiments serve to verify convergence, stability and accuracy, as well as a study of the performance of the proposed method on key features, such as mode conversion and Stoneley interface waves, arising from material interfaces. Numerical simulations for a curved interface and a curved cavity was also performed and the numerical results are in agreement with the theory. An open question in the SBP–SAT context is how to treat grid refinement for equations containing second derivatives in space. A study involving grid refinement for first order systems show good results using SBP preserving interpolation operators [18] but problems have been encountered when applied to higher order equations [12]. For the method to be more efficient grid refinement is a key feature. Also, an analytical study of the behavior of simulations involving Stoneley interface waves in materials with λ  μ is the topic of a coming study, as it was seen in [16] that the similar Rayleigh waves suffer increasingly from truncation errors in such materials. Acknowledgements The work of the first author was supported by King Abdullah University of Science and Technology (KAUST) through a joint KAUST Academic Excellence Alliance (AEA) grant with Stanford. The first author also acknowledges the support of Eric M. Dunham during this work.

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

61

Fig. 12. The numerical solutions at times t = 6.4, t = 7.0 and t = 7.6.

Appendix A. Kronecker product Definition 4 (Kronecker products). Let A be an m-by-n matrix and B be a p-by-q matrix. Then A ⊗ B, the Kronecker product of A and B, is the (mp)-by-(nq) matrix

62

K. Duru, K. Virta / Journal of Computational Physics 279 (2014) 37–62

⎛

⎞ . . . a1n B ... ... ... ⎟ ⎜ ... A⊗B =⎝ ⎠. ... ... ... ... am1 B am2 B . . . amn B a11 B

a12 B

The following properties hold for Kronecker products 1. Assume that the products AC and B D are well defined then ( A ⊗ B )(C ⊗ D ) = ( AC ) ⊗ ( B D ) 2. If A and B are invertible, then ( A ⊗ B )−1 = A −1 ⊗ B −1 3. ( A ⊗ B ) T = A T ⊗ B T References [1] L. Cagniard, E.A. Flinn, C.H. Dix, Reflection and Refraction of Progressive Seismic Waves, McGraw–Hill Book Company Inc., 1962. [2] M. Carpenter, D. Gottlieb, S. Abarbanel, Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: methodology and application to high-order compact schemes, J. Comput. Phys. 111 (1994). [3] M.H. Carpenter, J. Nordström, D. Gottlieb, A stable and conservative interface treatment of arbitrary spatial accuracy, J. Comput. Phys. 148 (1999) 341–365. [4] G. Cohen, S. Fauqueux, Mixed spectral finite elements for the linear elasticity system in unbounded domains, SIAM J. Sci. Comput. 26 (3) (2005). [5] M. Dehghannejad, A. Malehmir, C. Juhlin, P. Skyttä, 3D constraints and finite-difference modeling of massive sulfide deposits: the Kristineberg seismic lines revisited, northern Sweden, Geophysics 77 (2012) 69–79. [6] K. Duru, A perfectly matched layer for the time-dependent wave equation in heterogeneous and layered media, J. Comput. Phys. 257 (2014) 757–781. [7] K. Duru, G. Kreiss, K. Mattsson, Accurate and stable boundary treatments for elastic wave equations in second order formulation, SIAM J. Sci. Comput. (2014), in press. [8] K. Duru, K. Mattsson, G. Kreiss, Stable and conservative time propagators for second order hyperbolic systems, Tech. Rep., Div. Sc. Comp., Dept. of Infor. Tech., Uppsala University, 2011, pp. 201–208. [9] W.M. Ewing, W.S. Jardetzky, F. Press, Elastic Waves In Layered Media, McGraw–Hill Book Company Inc., 1957. [10] J.C. Gilbert, P. Joly, Higher Order Time Stepping for Second Order Hyperbolic Problems and Optimal CFL Conditions, Num. Analys and Sci. Comp. for PDEs and Their Challenging Applicat., vol. 16, Springer, 2008. [11] K.F. Graff, Wave Motion in Elastic Solids, Dover Publications, 1991. [12] M. Gustafsson, A. Nissen, K. Korrmann, Stable difference methods for block – structured adaptive grids, Technical report 2011-022, Department of Information Technology, Uppsala University, 2011. [13] J. Kozdon, E. Dunham, J. Nordström, Simulation of dynamic earthquake ruptures in complex geometries using high-order finite difference methods, J. Sci. Comput. (2012), http://dx.doi.org/10.1007/s10915-012-9624-5. [14] M. Käser, M. Dumbser, An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes. I. The two-dimensional isotropic case, Geophys. J. Int. 167 (2006). [15] H.-O. Kreiss, J. Oliger, Comparison of accurate methods for the integration of hyperbolic equations, Tellus 24 (1972). [16] H.-O. Kreiss, N.A. Petersson, Boundary estimates for the elastic wave equation in almost incompressible materials, SIAM J. Numer. Anal. 50 (3) (2012). [17] K. Mattsson, Summation by parts operators for finite difference approximations of second-derivatives with variable coefficients, J. Sci. Comput. 51 (2012) 650–682. [18] K. Mattsson, M. Carpenter, Stable and accurate interpolation operators for high-order multiblock finite difference methods, SIAM J. Sci. Comput. 32 (4) (2010). [19] K. Mattsson, F. Parisi, Stable and accurate second-order formulation of the shifted wave equation, Commun. Comput. Phys. 7 (2010) 103–137. [20] K. Mattsson, J. Nordström, Summation by parts operators for finite difference approximations of second derivatives, J. Comput. Phys. 199 (2004). [21] K. Mattsson, F. Ham, G. Iaccarino, Stable and accurate wave propagation in discontinuous media, J. Comput. Phys. 218 (2006) 8753–8767. [22] K. Mattsson, M. Svärd, M. Shoeybi, Stable and accurate schemes for the compressible Navier–Stokes equations, J. Comput. Phys. 227 (2008) 2293–2316. [23] N.A. Petersson, B. Sjögreen, Stable grid refinement and singular source discretization for seismic wave simulations, Commun. Comput. Phys. 8 (2010) 1074–1110. [24] F. Press, J. Oliver, M. Ewing, Seismic model study of refraction from a layer of finite thickness, Geophysics VXIX (3) (1954). [25] Rayleigh Lord, On waves propagated along the surface along the plane surface of an elastic solid, Proc. Lond. Math. Soc. 17 (1885). [26] R.A. Stephen, Solutions to range-dependent benchmark problems by the finite-difference method, J. Acoust. Soc. Am. 87 (1990). [27] R. Stoneley, Elastic waves at the surface of separation of two solids, Proc. R. Soc. Lond. A 6 (1924). [28] M. Svärd, J. Nordström, On the order of accuracy for difference approximations of initial-boundary value problems, J. Comput. Phys. 218 (October 2006) 333–352. [29] J. Virieux, P–SV wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics 51 (1986) 889–901. [30] K. Virta, G. Kreiss, Surface waves in almost incompressible elastic materials, Technical report, arXiv:1309.3863 [physics.geo-ph]. [31] K. Virta, K. Mattsson, Acoustic wave propagation in complicated geometries and heterogeneous media, J. Sci. Comput. (2014), http://dx.doi.org/ 10.1007/s10915-014-9817-1. [32] V. Lisitsa, D. Vishnevskiy, Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity, Geophys. Prospect. 58 (2010) 619–635. [33] E.H. Saenger, N. Gold, S.A. Shapiro, Modeling the propagation of the elastic waves using a modified finite-difference grid, Wave Motion 31 (2000) 77–92. [34] A.R. Levander, Fourth-order finite-difference P-SV seismograms, Geophysics 53 (1988) 1425–1436. [35] J. Virieux, P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics 51 (1986) 889–901. [36] B. Sjögreen, N.A. Petersson, A fourth order accurate finite difference scheme for the elastic wave equation in second order formulation, J. Sci. Comput. 52 (1) (2012) 17–48. [37] S. Nilsson, N.A. Petersson, B. Sjögreen, H.-O. Kreiss, Stable difference approximations for the elastic wave equation in second order formulation, SIAM J. Numer. Anal. 45 (5) (2007) 1902–1936. [38] D. Appelö, N.A. Petersson, A stable finite difference method for the elastic wave equation on complex geometries with free surfaces, Commun. Comput. Phys. 5 (2009) 84–107.