MA471 Fall 2003

1
Lecture 22

• MA471 Fall 2003

2
Advection Equation

• Recall the 2D advection equation
• We will use a Runge-Kutta time integrator and
  spectral representation in space.

3
Periodic Data

• Lets assume we are given N values of a function
  f at N data points on the unit interval.
• For instance N8 on a unit interval

0 1/8 2/8 3/8 4/8 5/8
6/8 7/8 8/8
4
Discrete Fourier Transform

• We can seek a trigonometric interpolation of a
  function f as a linear combination of N (even)
  trigonometric functions
• Such that

5
Transform

• The interpolation formula defines a linear system
  for the unknown fhat coefficients

Or
6
Code for the DFT
7
Code for Inverse DFT
8
Fast Fourier Transform

• See handout

9
Spectral Derivative

• We can differentiate the interpolant by

10
Detail

• We note that the derivative of the k(N/2) mode
• is technically
• However, as we show on the next slide this mode
  has turning points at all the data points.

11
Real Component of N/2 Mode
i.e. the slope of the k(N/2) mode is zeroat all
the 8 points..
12
Modified Derivative Formula

• So we can create an alternative symmetric
  derivative formula

13
Spectral Differentiation Scheme

1. Use an fft to compute
2. Differentiate in spectral space. i.e. compute

14
cont

• 3) Then
• 4) Summary
• a) fft transform data to compute coefficients
• b) scale Fourier coefficients
• c) inverse fft (ifft) scaled coefficients

15
Final Twist

• Matlab stores the coefficients from the fast
  Fourier transform in a slightly odd order
• The derivative matrix will now be a matrix with
  diagonal entries

16
Spectral Differentiation Code
Typo See corrected Code on webpage

1. DFT data
2. Scale Fourier coefficients
3. IFT scaled coefficients

17
Two-Dimensional Fourier Transform

• We can now construct a Fourier expansion in two
  variables for a function of two variables..
• The 2D inverse discrete transform and discrete
  transform are

18
Advection Equation

• Recall the 2D advection equation
• We will use a Runge-Kutta time integrator and
  spectral representation in space.

19
Runge-Kutta Time Integrator

• We will now discuss a particularly simple
  Runge-Kutta time integrator introduced by
  Jameson-Schmidt-Turkel
• The idea is each time step is divided into s
  substeps, which taken together approximate the
  update to sth order.

20
Side note Jameson-Schmidt-Turkel Runge-Kutta
Integrator

• Taylors theorem tells us that
• We will compute an approximate update as

21
JST Runge-Kutta

• The numerical scheme will look like
• We then factorize the polynomial derivative term

22
Factorized Scheme
23
JST Advection Equation

• We now use the PDE definition

24
Now Use Spectral Representation

• A time step now consists of s substages.
• Each stage involves
• Fourier transforming the ctildeusing a fast
  Fourier transform (fft)
• Scaling the coefficients to differentiate in
  Fourier space
• Transforming the derivatives back tophysical
  values at the nodes by inverse fast Fourier
  transform (ifft).
• Finally updating ctilde according tothe
  advection equation.
• At the end we update the concentration.

25
Matlab Implementation

• First set up the nodes and the differentiation
  scalings.

26

• 24) Compute dt
• 26) Initiate concentration
• 28-29) compute advection vector
• 32-45) full Runge-Kutta time step
• 35) fft ctilde
• 37) x- and y-differentiation in Fourier space
• 40-41) transform back to physical values of
  derivatives
• 43) Update ctilde

27
(No Transcript)
28
FFT Resources

• Check out http//www.fftw.org/ for the fastest
  Fourier transform in the West
• Great list of links
• http//www.fftw.org/links.html
• FFT for irregularly spaced data
• http//www.math.mu-luebeck.de/potts/nfft/

29
Lab Exercise

1. Download example DFT code from website.
2. Benchmark codes as a function of N
3. Write your own version using an efficient FFT
   (say fftw) which
4. FFTs data
5. Scales coefficients
6. IFFTs scaled coeffients
7. Write an advection code using the JST Runge-Kutta
   time integrator.

This is the scalar version of Project 4
30
Upshot 2D Advection/DFTUsing MPI_Alltoall
Note using fft the compute time will drop by
approx. A factor of 30 (for N256)
31
Upshot 2D Advection/DFTUsing Non-Blocked Sends
Note I isend and irecv, then compute the
x-derivatives, then waitall, then compute
the y-derivatives, then isend/recv and waitall
(not the best option perhaps).
32
With FFT

• Replacing the DFT calls with FFT calls will
  radically reduce the computation time.
• It will be much harder to hide the communication
  time behind the compute time.
• Good luck ?