Transcript igel

Session:
Computational Wave
Propagation: Basic Theory
Igel H., Fichtner A., Käser M., Virieux J., Seriani G., Capdeville Y., Moczo P.
 The finite-difference method
 The pseudo-spectral method
The problem
The wave equation
 ( x) u(x,t)   x E(x) xu(x,t)
2
t
Displacement u, elastic parameters E(x) vary only in one direction.
Density  and Lame parameters l and m.
E ( x)  m ( x)
E ( x)  l ( x)  2 m ( x)
shear waves
P waves
with the corresponding velocities
m
vS 

shear waves
vP 
l  2m

P waves
… we disregard …
 3-dimensionality
 External forces
 Viscoelastic behaviour
 Anisotropic behaviour
 Non-linear elasticity
 Etc.
… and even further simplify and re-write to …
… motion of a string …
 ( x) u(x,t)   x m (x) xu(x,t)
2
t
Displacement u, shear modulus m(x), density 
vary only in one direction x.
The velocity-stress formulation
1
 t u ( x, t ) 
 x ( x, t )
 ( x)
t ( x, t )  m ( x) xu( x, t )
Extrapolation
usually low-order
stress-strain
Generalized Interpolation
Low-order to spectral accuracy
Other forms of the elastic wave equation
 Strong form (see Moczo et al.)




Displacement – stress
Displacement – velocity – stress
Velocity – stress
Displacement
 Weak form (more later)
 Projection on test functions
Numerical methods






Finite differences
Pseudospectral methods
Finite (spectral) elements (Fichtner)
Discontinuous Galerkin (Käser)
Frequency domain methods (Virieux)
Cell methods (Seriani)
(Boundary integral methods, lattice solids, discrete
wavenumber, finite volumes, etc).
Finite Differences
f ( x  dx )  f ( x)
x f 
dx
forward difference
f ( x)  f ( x  dx )
x f 
dx
backward difference
f ( x  dx )  f ( x  dx )
x f 
2dx
centered difference


c
... or even …
f ( x  dx / 2)  f ( x  dx / 2)
x f 
dx
c
… in space and for time …
f (t  dt / 2)  f (t  dt / 2)
t f 
dt
c
This turns out to be a useful approximation!
Discrete space-time
Space and time
discretization
Regular, 1D
t  ldt
x  mdx
t  dt  (l  1)dt
1
t  dt / 2  (l  )dt
2
x  dx  ( m  1)dx
1
x  dx / 2  ( m  )dx
2
Our first FD wave algorithm
t  ldt
x  mdx
Space discretization
Regular, 1D
 u

1 

dt
m
dx
l
l 1
l 1/ 2
l 1/ 2
 m1/ 2   m1/ 2
um1  um
 mm
dt
dx
u
l 1/ 2
m
l 1/ 2
m
l
m1/ 2
l
m1/ 2
… and the explicit extrapolation scheme …
Source term
Do for all times …


dt l
dt l
l
u
 u 
 m1/ 2   m1/ 2 
fm
dx m
m
dt l 1/ 2 l 1/ 2
l
l 1
 m1/ 2   m1/ 2  mm um1  um 
dx
l 1/ 2
m
Stop!
l 1/ 2
m
… staggered grids … ?
the only unknown

Time
l+1/2
l
l-1/2




m-1/2 m m+1/2
Space
Results
Velocity 5 km/s
3000
2500
Time (s)
2000
1500
1000
500
0
0
1000
2000
3000
4000
Distance (km)
5000
6000
Amplitude
Numerical dispersion
Time
Von Neumann Analysis
plane wave analysis
continuous
f ( x, t )  A exp(ikx  it )
flm  A exp(ikmdx  ildt )
Conditional stability
 dt 
vP,S 
 1
 dx 
discrete
Numerical dispersion
phase velocity - group velocity
Theoretical prediction of phase
and group velocities as a function
of wavelength.
In 2-D, 3-D this effect depends on
direction (-> grid anisotropy)
Well, let us have a closer look at the space
derivatives (because there is not so much we
can do about the time derivative …)


ikx
 x f ( x)   x   F (k )e dk 
 


   ikF (k )e

ikx
dk
The pseudospectral approach

 x f ( x)    ikF (k )e
ikx
dk

in the discrete world that means, we can …
DFT f m   Fn
Discrete Spectrum
Fn  (ik ) Fn
IDFT(ik ) Fn    x f m
Exact derivative
PS Algorithm (displacement)
1
 x μ(x) xu(x,t)  2u(x,t) u(x,t dt)
u(x,t dt)  2
dt ρ(x)
Derivatives calculated by DFT (forward and
inverse transform and multiplication with –ik)
…. cool algorithm, but …
… some more thoughts on the spectral derivative
operator „-ik“ … using the convolution theorem …

 x f ( x)    ikF (k )e

ikx

dk   D(k ) F (k )e


  d ( x' ) f ( x  x' )dx

Differential operator
in space domain
The function we
want to differentiate
ikx
dk
„Convolution“
f
d
0
0
1
0
-2
1
1
0
0
1
-2
1
0
1
0
0
1
-2
1
1
0
1
-2
1
0
-2
0
-2
0
1
0
1
0
1
1
d*f
1
1
0
0
0
0
1
0
0
1
-2
1
0
The numerical wave number as a
function of true wave number (ik as
a function of k)
This is the derivative operater d(x)
and it s discrete version (dots)
Finite difference - Pseudospectral
2-4-6 point
Taylor
coefficients
FD
Length of (convolutional) operator
Increasing accuracy
Increasing floating point operations
nx-point
Fourier
(Chebysheff)
coefficients
PS
FD in real 3-D applications
 Cartesian Geometry
 Spherical coordinates
 Sections
 Axisymmetric
 Rheologies
 Viscolastic
 Anisotropic
 Poroelastic
 Dynamic rupture
Finite Differences - Pros and Cons
PROS
CONS


 Boundary conditions difficult to
implement
 Requires large number of grid
points per wavelength (in
particular surface waves)
 Large memory requirement
 Inefficient for models with
strong velocity variations
 Topography not easily
implemented



simple theory
Explicit scheme – no matrix
inversion
Easy parallelization
Easy model generation (on
regular grids)
Easy adaptation to specific
problems
Pseudospectral Method - Pros and Cons
PROS
CONS



 Global communication scheme
(inefficient parallelization)
 Irregular grid points
(Chebyshev) –> stability
problems
 Boundary conditions difficult
for Fourier Method



Beautiful!
Exact space derivatives
Explicit scheme – no matrix
inversion
Centred scheme (anisotropy)
Accurate implementation of
boundary condition
(Chebyshev)
Memory efficient (less points
per wavelength w.r.t. FD)
Final thoughts


Some of the fundamental concepts of computational
seismology can be understood from FD schemes (e.g.,
stability, dispersion)
Each QUEST researcher should be able to




Understand what time steps and grid distances are
appropriate for a specific Earth velocity models
Understand what affects the accuracy of wave simulations
(np/l; operators used, wave type, complexity of velocity
model, propagation distance, etc.)
Know the possible traps of (community) algorithms (black
boxes)
Know how to benchmark simulation codes (for accuracy)