320 likes | 568 Vues
Lecture 22. MA471 Fall 2003. Advection Equation. Recall the 2D advection equation: We will use a Runge-Kutta time integrator and spectral representation in space. Periodic Data. Let’s assume we are given N values of a function f at N data points on the unit interval.
E N D
Lecture 22 MA471 Fall 2003
Advection Equation • Recall the 2D advection equation: • We will use a Runge-Kutta time integrator and spectral representation in space.
Periodic Data • Let’s assume we are given N values of a function f at N data points on the unit interval. • For instance N=8 on a unit interval: 0 1/8 2/8 3/8 4/8 5/8 6/8 7/8 8/8
Discrete Fourier Transform • We can seek a trigonometric interpolation of a function f as a linear combination of N (even) trigonometric functions: • Such that:
Transform • The interpolation formula defines a linear system for the unknown fhat coefficients: Or:
Fast Fourier Transform • See handout
Spectral Derivative • We can differentiate the interpolant by:
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.
Real Component of N/2 Mode i.e. the slope of the k=(N/2) mode is zeroat all the 8 points..
Modified Derivative Formula • So we can create an alternative symmetric derivative formula:
Spectral Differentiation Scheme • Use an fft to compute: • Differentiate in spectral space. i.e. compute:
cont 3) Then: 4) Summary: a) fft transform data to compute coefficients b) scale Fourier coefficients c) inverse fft (ifft) scaled coefficients
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:
Spectral Differentiation Code Typo: See corrected Code on webpage • DFT data • Scale Fourier coefficients • IFT scaled coefficients
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:
Advection Equation • Recall the 2D advection equation: • We will use a Runge-Kutta time integrator and spectral representation in space.
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 s’th order.
Side note: Jameson-Schmidt-Turkel Runge-Kutta Integrator • Taylor’s theorem tell’s us that • We will compute an approximate update as:
JST Runge-Kutta • The numerical scheme will look like • We then factorize the polynomial derivative term:
JST + Advection Equation • We now use the PDE definition +
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.
Matlab Implementation • First set up the nodes and the differentiation scalings.
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
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/
Lab Exercise • Download example DFT code from website. • Benchmark codes as a function of N • Write your own version using an efficient FFT (say fftw) which: • FFT’s data • Scales coefficients • IFFT’s scaled coeffients • Write an advection code using the JST Runge-Kutta time integrator. This is the scalar version of Project 4
Upshot: 2D Advection/DFTUsing MPI_Alltoall Note: using fft the compute time will drop by approx. A factor of 30 (for N=256)
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).
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