Transcript Pseudospectral Methods - uni

Pseudospectral Methods





What is a pseudo-spectral Method?
Fourier Derivatives
The Acoustic Wave Equation with the Fourier Method
Comparison with the Finite-Difference Method
Dispersion and Stability of Fourier Solutions
The goal of this lecture is to shed light at one end of the axis of FD
(or convolutional) type differential operators. When one uses all
information of a space-dependent field to estimate the derivative at a
point one obtains spectral accuracy.
Pseudospectral methods
1
What is pseudospectral?
Spectral solutions to time-dependent PDEs are formulated
in the frequency-wavenumber domain and solutions are
obtained in terms of spectra (e.g. seismograms). This
technique is particularly interesting for geometries where
partial solutions in the -k domain can be obtained
analytically (e.g. for layered models).
In the pseudo-spectral approach - in a finite-difference like
manner - the PDEs are solved pointwise in physical space
(x-t). However, the space derivatives are calculated using
orthogonal functions (e.g. Fourier Integrals, Chebyshev
polynomials). They are either evaluated using matrixmatrix multiplications, fast Fourier transform (FFT), or
convolutions.
Pseudospectral methods
Fourier Derivatives
.. let us recall the definition of the derivative using Fourier integrals ...


ikx
 x f ( x)   x   F (k )e dk 
 


   ikF (k )e ikxdk

... we could either ...
1) perform this calculation in the space domain by convolution
2) actually transform the function f(x) in the k-domain and back
Pseudospectral methods
The acoustic wave equation
let us take the acoustic wave equation with variable density
1

1
2

t p   x   x p 
2

c


the left hand side will be expressed with our
standard centered finite-difference approach
1

1

 p(t  dt)  2 p(t )  p(t  dt)   x   x p 
2
2
c dt


... leading to the extrapolation scheme ...
Pseudospectral methods
The Fourier Method: acoustic wave propagation
1


p(t  dt)  c dt  x   x p 
 2 p(t )  p(t  dt)



2
2
where the space derivatives will be calculated using the Fourier Method.
The highlighted term will be calculated as follows:
ˆ n  ik P
ˆ n  FFT1   P n
Pjn  FFT  P

 
x j
multiply by 1/
n
n
 1 ˆ
 1 ˆ
1
1
n
1
n
 x Pj  FFT    x P   ik   x P   FFT   x   x Pj 







... then extrapolate ...
Pseudospectral methods
… and in 3D …
p (t  dt) 
 1

1

1


c dt   x   x p    y   y p    z   z p  





 
 2 p (t )  p (t  dt)
2
2
.. where the following algorithm applies to each space dimension ...
ˆ n  ik P
ˆ n  FFT1   P n
Pjn  FFT  P

 
x j
n
n
 1 ˆ
 1 ˆ
1
1
n
1
n
 x Pj  FFT    x P   ik   x P   FFT   x   x Pj 







Pseudospectral methods
… FD in Matlab …
let us compare the core of the algorithm - the calculation of the derivative
(Matlab code)
function
df=fder1
%
fDER1D(f,dx,nop
%
second
derivati
nx=max(size(f));
n2=(nop-1)/2;
if
if
nop==3;
nop==5;
d=[1
d=[-1/
df=[1:nx]*0;
for
i=1:nop;
df=df+d(i).*cshif
end
Pseudospectral methods
… and as PS …
... and the first derivative using FFTs ...
function
df=
%
SDER1D(f,d
nx=max(size(
%
initialize
kmax=pi/dx;
dk=kmax/(nx/
for
i=1:nx/2
k=sqrt(-1)*k
%
FFT
and
IF
ff=fft(f);
f
.. simple and elegant ...
Pseudospectral methods
Dispersion and Stability
i ( kjdxndt )
p e
n
j
2 i ( kjdxndt )
 p  k e
2
x
n
j
4
2 dt i ( kjdxndt )
 p   2 sin
e
dt
2
2
t
n
j
2
dt
k
sin
cdt
2
Pseudospectral methods
Dispersion and Stability
2
dt
k
sin
cdt
2
2
1 kcdt
  sin (
)
dt
2
What are the consequences?
a) when dt << 1, sin-1 (kcdt/2) kcdt/2 and w/k=c
-> practically no dispersion
b) the argument of asin must be smaller than one.
kmaxcdt
1
2
cdt / dx  2 /   0.636
Pseudospectral methods
Results @ 10Hz
Source time function
1
1
0.5
0.5
0
-0.5
200
400
600
0
Example of acoustic 1D wave simulation.
3 point - 2 order; T = 6.6 s, E
FD3 -point operator
red-analytic;
blue-numerical; green-difference
1
Pseudospectral methods
0
0
Results @ 10Hz
Source time function
1
1
0.5
0.5
0
-0.5
Pseudospectral methods
0
0
200
400
600
Example of acoustic 1D wave simulation.
point - 2 order; T = 7.8 s,
FD 5 -point 5
operator
red-analytic; blue-numerical; green-difference
1
Results @ 10Hz
Source time function
1
1
0.5
0.5
0
-0.5
200
400
600
0
- 2 order; T = 35 s, Er
Example of acousticFourier
1D wave simulation.
Fourier operator
1
red-analytic; blue-numerical; green-difference
Pseudospectral methods
0
0
Results @ 20Hz
Source time function
1
1
0.5
0.5
0
-0.5
200
400
600
0
- 2 order; T = 7.8 s, E
Example of acoustic3
1Dpoint
wave simulation.
FD3 -point operator
1
red-analytic; blue-numerical; green-difference
Pseudospectral methods
0
0
Results @ 20Hz
Source time function
1
1
0.5
0.5
0
-0.5
200
400
600
Example of acoustic 1D wave simulation.
5 point - 2 order; T = 7.8 s,
FD 5 -point operator
red-analytic;
blue-numerical; green-difference
1
Pseudospectral methods
0
0
Results @ 20Hz
Source time function
1
1
0.5
0.5
0
-0.5
200
400
600
0
Example of acoustic
1D wave simulation.
Fourier
- 2 order; T = 34 s, Erro
Fourier operator
1
red-analytic;
blue-numerical; green-difference
Pseudospectral methods
0
0
Computational Speed
Difference (%) between numerical and analytical solution
as a function of propagating frequency
160
140
120
100
3 point
5 point
Fourier
80
60
40
20
0
5 Hz
Pseudospectral methods
10 Hz
20 Hz
Simulation time
5.4s
7.8s
33.0s
Numerics and Green‘s Functions
The concept of Green’s Functions (impulse responses) plays an
important role in the solution of partial differential equations. It is also
useful for numerical solutions. Let us recall the acoustic wave equation
t2 p  c 2 p
with  being the Laplace operator. We now introduce a delta source in
space and time
t2 p   ( x) (t )  c 2 p
the formal solution to this equation is
1  (t  x / c)
p( x, t ) 
4c 2
x
(Full proof given in Aki and Richards, Quantitative Seismology, Freeman+Co, 1981, p. 65)
Pseudospectral methods
Numerical Green‘s functions
Fourier Method
10
10
10
9
9
9
8
8
8
7
7
7
6
6
6
5
5
5
4
4
4
3
3
3
2
2
2
1
1
1
0
500
Pseudospectral methods
1000
1500
0
500
1000
0
1500 500
1000
1500
Impulse response (numerical convolved with source
5 point operator
Impulse response (analytical) concolved with source
Frequency increases
3 point operator
Pseudospectral Methods - Summary
The Fourier Method can be considered as the limit of the finite-difference
method as the length of the operator tends to the number of points
along a particular dimension.
The space derivatives are calculated in the wavenumber domain by
multiplication of the spectrum with ik. The inverse Fourier transform
results in an exact space derivative up to the Nyquist frequency.
The use of Fourier transform imposes some constraints on the
smoothness of the functions to be differentiated. Discontinuities lead
to Gibb’s phenomenon.
As the Fourier transform requires periodicity this technique is particular
useful where the physical problems are periodical (e.g. angular
derivatives in cylindrical problems).
Pseudospectral methods play a minor role in seismology today but the
principal of spectral accuracy plays an important role in spectral
element methods
Pseudospectral methods