Lecture 22Advection EquationPeriodic DataDiscrete Fourier TransformTransformCode for the DFTCode for Inverse DFTFast Fourier TransformSpectral DerivativeDetailReal Component of N/2 ModeModified Derivative FormulaSpectral Differentiation SchemecontFinal TwistSpectral Differentiation CodeTwo-Dimensional Fourier TransformAdvection EquationRunge-Kutta Time IntegratorSide note: Jameson-Schmidt-Turkel Runge-Kutta IntegratorJST Runge-KuttaFactorized SchemeJST + Advection EquationNow Use Spectral RepresentationMatlab ImplementationFFT ResourcesLab ExerciseUpshot: 2D Advection/DFTUsing MPI_AlltoallUpshot: 2D Advection/DFTUsing Non-Blocked SendsWith FFTLecture 22MA471 Fall 2003Advection Equation• Recall the 2D advection equation:• We will use a Runge-Kutta time integrator and spectral representation in space.0xyCCCaatxy∂∂∂++=∂∂∂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/8Discrete Fourier Transform• We can seek a trigonometric interpolation of a function f as a linear combination of N (even) trigonometric functions:• Such that:()()/22/2 1ˆkNik xkkNIf x f eπ==− +=∑() ()()/22/2 1ˆfor 1,...,where jkNik xjj kkNjfx Ifx fe j NjxNπ==− +== ==∑Transform• The interpolation formula defines a linear system for the unknown fhat coefficients:()/22/2 121ˆ for 1,...,1ˆ for k 1,...,22jjkNik xjkkNjNik xkjjffejNNNffeNππ==− +=−=====−+∑∑Or:()()ˆˆifftfft==ffffCode for the DFTCode for Inverse DFTFast Fourier Transform• See handoutSpectral Derivative• We can differentiate the interpolant by:()()/22/2 1ˆkNikxkkNIf x f eπ==− +=∑() ( )()/22/2 1ˆ2kNik xkkNdIfxikfedxππ==− +=∑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.222ˆiNxNefπ222ˆ22iNxNiNefππReal Component of N/2 Modei.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:() ()()()/2 12/2 11ˆ22kNik xkkNdIfxikfedxπππ=−=− +≅∑Spectral Differentiation Scheme1) Use an fft to compute:2) Differentiate in spectral space. i.e. compute:ˆ for 1, 2,..., 1,0,1,...,22 2kNN Nfk =− + − + − ()ˆ2 for 1,..., 122ˆ0 for 2kkNNikf kgNkπ =−+ − ==cont3) Then:4) Summary:a) fft transform data to compute coefficientsb) scale Fourier coefficientsc) inverse fft (ifft) scaled coefficients()()()/2/2 1ˆkNikxkkNdIfxgedx==− +≅∑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:01 11122222ˆˆ ˆ ˆ ˆ ˆ, ,.., , , , ,...,ˆNNNNff fffff−−−+−+2 0,1,.., 1, , 1, 2,..., 1222NNNiπ−−+−+ −0Spectral Differentiation CodeTypo: See correctedCode on webpage1) DFT data2) Scale Fourier coefficients3) IFT scaled coefficientsTwo-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:()()/2/222/2 1 /2 122211ˆ for , 1,...,1ˆ for j,k 1,...,22nmnmjNkNij x ik ynm jkkN jNnNmNij x ik yjk nmnmffeenmNNNffeeNππππ===− + =− +==−−======−+∑∑∑∑()()()/2/222/2 1 /2 1ˆ,jNkNij x ik yjkkN jNIf x y f e eππ===− + =− +=∑∑Advection Equation• Recall the 2D advection equation:• We will use a Runge-Kutta time integrator and spectral representation in space.0xyCCCaatxy∂∂∂++=∂∂∂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:() ()()()[]2211**11 ...1! 2! ! for some ,1!sssssdt d dt d dt dCt dt Ctdt dt s dtdt d Cttttdtsdt+++ +=+ + ++ +∈++() ()221 ...1! 2! !ssdt d dt d dt dCt dt Ctdt dt s dt +≅+ + ++ JST Runge-Kutta• The numerical scheme will look like• We then factorize the polynomial derivative term:() ()222211 ...1! 2! !1 ...1! 2! !ssssn ndt d dt d dt dCt dt Ctdt dt s dtdt d dt d dt dCCdt dt s dt+ +≅+ + ++ ⇒=+ + ++ Factorized Scheme2211 ...1! 2! !ssn ndt d dt d dt dCCdt dt s dt+ =+ + ++ n+1..121nnn n ndt d dt d dt d dt dCCC C C Csdt s dt s dt dt=+ + + +−−nnn+1Set C=Cfor : 1:1 C CendC=Cksdt dCkdt=−←+JST + Advection Equation• We now use the PDE definition nnn+1Set C=Cfor : 1:1 C CendC=Cksdt dCkdt=−←+0xyCCCaatxy∂∂∂++=∂∂∂+nnn+1Set C=Cfor : 1:1 C CendC=Cxyksdt C Caakx y=−∂∂←− +∂∂Now Use Spectral Representation()()()()jjSet for : 1:1ˆ = fftˆˆ ˆˆ ˆ ifftˆ ifft for j=1,..,Nend=xxyyxxyyxx yyjj jjksdtk==−====←− +ccccdDcdDcddddcc adadccA time step now consists of s substages.Each stage involves:a) Fourier transforming the ctildeusing a fast Fourier transform (fft)b) Scaling the coefficients to differentiate in Fourier spacec) Transforming the derivatives back tophysical values at the nodes by inverse fast Fourier transform (ifft).d) Finally updating ctilde according tothe advection equation.At the end we update the concentration.MatlabImplementation• First set up the nodes and the differentiation scalings.24) Compute dt26) Initiate concentration28-29) compute advection vector32-45) full Runge-Kutta time step35) fft ctilde37) x- and y-differentiation in Fourier
View Full Document