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 m1/ 2 m1/ 2 um1 um mm dt dx u l 1/ 2 m l 1/ 2 m l m1/ 2 l m1/ 2 … and the explicit extrapolation scheme … Source term Do for all times … dt l dt l l u u m1/ 2 m1/ 2 fm dx m m dt l 1/ 2 l 1/ 2 l l 1 m1/ 2 m1/ 2 mm um1 um 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 it ) flm A exp(ikmdx ildt ) 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)