A fourth order accurate finite difference scheme for the elastic - PowerPoint PPT Presentation

A fourth order accurate finite difference scheme for the elastic wave equation in second order formulation 1 B. Sj¨ ogreen and N.A.Petersson Lawrence Livermore National Laboratory October 18, 2011 1Work performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-PRES-505131.

Computational seismology ρ u tt = div σ + f ◮ u = u ( x , y , z , t ) displacement vector ( u = ( u v w )). ◮ f = f ( x , y , z , t ) forcing = earthquake model ◮ stress tensor:   (2 µ + λ ) u x µ ( u y + v x ) µ ( u z + w x ) σ = µ ( u y + v x ) (2 µ + λ ) v y µ ( v z + w y )   µ ( u z + w x ) µ ( v z + w y ) (2 µ + λ ) w z ◮ ρ = ρ ( x , y , z ), µ = µ ( x , y , z ), λ = λ ( x , y , z ) mtrl. prop.

Domain and wave types Computational domain ◮ Traction free boundary condition at surface. � ◮ Pressure wave with speed c p = (2 µ + λ ) /ρ . ◮ Shear wave with speed c s = � µ/ρ . √ ◮ Wave speed ratio c p / c s > 2. ◮ Rayleigh waves on surface, slower than P- and S-waves.

Topography handled by curvilinear grid Grid refinement for depth varying wave speeds.

Resolution requirements h = min c s Pf ◮ Grid spacing h ◮ Points per shortest wavelength P ◮ Highest frequency f ◮ Material shear wave speed c s Typical values: f = 10 Hz, c s = 300 m/s, P = 15 (second order), P = 7 (fourth order), gives h = 2 m (2nd order) h = 4 . 29 m (4th order) Domain size 200 km → 100,000 pts/dimension (2nd) 46,620 (4th)

Objective This work (new): 4th order accurate energy conserving method. Previous work: 2nd order accurate energy conserving method. Extensions to ◮ Curvilinear grids ◮ Far field boundaries ◮ Mesh refinement ◮ Viscoelastic model

Energy conserving methods for the elastic wave equation E n discrete energy at t n , integral over space, conserved when f = 0 E n = E n − 1 = . . . = E 0 . Compatibility with norm, c 1 || u n || h ≤ E n ≤ c 2 || u n || h gives stability, || u n || h ≤ E n / c 1 = . . . = E 0 / c 1 ≤ c 2 / c 1 || u 0 || h ◮ Stability for inhomogeneous material, real b.c., any c p / c s . ◮ Stable for long time integration ◮ Dissipation free ◮ Robust code, no numerical parameters to tune, but must be careful to not introduce unresolved frequencies

Energy estimate gives long time stability Standard stability gives convergence on 0 < t < T with T fixed.

Example: Wave equation with mixed derivative term ( x , y ) ∈ [0 , 1] 2 , t > 0 u tt = (2 au x + au y ) x + ( au x + 2 au y ) y , a = a ( x , y ) > 0 variable coefficient. Boundary conditions: u = 0 at x = 0 2 u x + u y = 0 at x = 1 u ( x , y , t ) = u ( x , y + 1 , t ) (periodic in y )

Energy estimate 1 d || u t || 2 + ( u x , au x ) + ( u x + u y , a ( u x + u y )) + ( u y , au y ) � � = 0 2 dt (Note: Non-negative terms give L 2 estimate) Derived by partial integration: 1 d dt || u t || 2 = ( u t , u tt ) = . . . = 2 − 1 d dt (( u x , 2 au x ) + ( u x , au y ) + ( u y , au x ) + ( u y , 2 au y )) + B . T 2 Energy terms: ( u x , au x ) + ( u x + u y , a ( u x + u y )) + ( u y , au y ) B . T . = u t a (2 u x + u y ) | x =1 − u t a (2 u x + u y ) | x =0 zero by b.c.

Discretization Cartesian grid with constant spacing h . Centered finite difference operators ∂ u ( x i ) /∂ x → D 0 u i , i = 1 , . . . , N satisfying summation-by-parts ( u , D 0 v ) h = − ( D 0 u , v ) h + u N v N − u 1 v 1 in a discrete, weighted, scalar product ( u , v ) h . Further notation: D + u i = ( u i +1 − u i ) / h , D − u i = ( u i − u i − 1 ) / h . In two dimensions: D ( x ) 0 u i , j and D ( y ) 0 u i , j .

Discretization ( au y ) x ≈ D ( x ) 0 ( aD ( y ) 0 u ) and ( au x ) x ≈ D ( x ) 0 ( aD ( x ) 0 u ) same energy estimate as for PDE possible, but ◮ Energy not positive definite, norm estimate not possible. ◮ Boundary condition 2 u x + u y = 0 implicit. Second order method ( au x ) x ≈ D + ( a j − 1 / 2 D − u j ), where Energy estimate based on D + ( a j − 1 / 2 D − u j ) = D 0 ( a j D 0 u j ) − h 2 4 D + D − ( a j D + D − u j ) , Square completion with x - y terms Keeps energy pos. def. Use of ghost points, gives explicit discrete b.c. with no boundary modification of D + D − .

Fourth order accurate operator ( au x ) x ≈ G ( a , u ) j = D 0 ( a j D 0 u j )+ h 4 18 D + D − D + ( a j − 1 / 2 D − D + D − u j ) − h 6 144( D + D − ) 2 ( a j ( D + D − ) 2 u j ) + boundary modifications ◮ G is five point wide operator away from the boundary. ◮ D 0 SBP operator of order 4/2, needed for xy -derivatives. ◮ G also order 4/2. Boundary modified at j = 1 , . . . , 6. ◮ B.T.=0 in SBP is 4th order accurate b.c. → 4th order error. ◮ Boundary modification of ( D + D − ) 3 gives first order errors that can be made to cancel first order errors of D 0 ( aD 0 u ). ◮ Can expand G ( a , u ) j = � 8 � 8 k =1 β j , k , m a k u m , j = 1 , . . . , 6. m =1 Coefficient tensor β with 129 non-zero elements out of 384. ◮ G uses ghost points, D 0 does not.

4th order P-C time discretization gives energy conservation Can prove time discrete energy conservation: E n +1 / 2 = E n − 1 / 2 . Method stable (energy positive) for CFL < 1 . 3. No stiffness for high order.

Numerical examples Elastic wave equation, 2D ρ u tt = ((2 µ + λ ) u x ) x + ( λ v y ) x + ( µ v x ) y + ( µ u y ) y ρ v tt = ( µ v x ) x + ( µ u y ) x + ( λ u x ) y + ((2 µ + λ ) v y ) y 0 < x < L x , 0 < y < L y , t > 0. Initial data: u ( x , y , 0) and u t ( x , y , 0) given. Boundary data: y -periodic, with Dirichlet b.c. on x = L x and (2 µ + λ ) u x + λ v y = 0 x = 0 µ ( v x + u y ) = 0 x = 0

Energy test with random material λ ( x , y ) = 2( r 2 − 2) + θ ρ ( x , y ) = 4 + θ µ ( x , y ) = 2 + θ Random variable θ ∈ [0 , 1]. Approximate wave speed ratio r = c p / c s . Initial data also random numbers. Energy change per time step. Total > 220 , 000 steps. c p / c s arbitrarily large.

Rayleigh waves Surface waves at x = 0, solutions u s traveling wave in y and decaying as e − ax into the domain. µ , λ , and ρ constant.

34 seconds vs. 54 hours CPU time Error vs. CPU time µ = 0 . 001, error 10 − 4 need 34 seconds with 4th order scheme, 54 hours with 2nd order scheme.

Rayleigh waves

Summary and future directions ◮ 4th order accurate non-dissipative difference scheme, L 2 norm stable with heterogeneous material and boundary conditions. ◮ 4th order in both space and time. ◮ Significant savings in computational resources. ◮ High order second derivative approximation of ( µ ( x ) u x ) x , with norm stable boundary closure, useful in other applications. ◮ To be implemented into the 3D WPP solver. ◮ To be used in new solver for source and material inversion, using adjoint wave propagation. Reference [1] B. Sj¨ ogreen and N.A.Petersson, A fourth order finite difference scheme for the elastic wave equation in second order formulation , Lawrence Livermore National Laboratory, LLNL-JRNL-483427, (to appear in J.Scient.Comput.).