tag:blogger.com,1999:blog-5822805028291837738.post1012926326293474096..comments2020-07-05T05:46:40.545-04:00Comments on Various Consequences: FFTW Discrete Cosine Transform DerivativeJoshua Stultshttp://www.blogger.com/profile/03506970399027046387noreply@blogger.comBlogger19125tag:blogger.com,1999:blog-5822805028291837738.post-53820978692561190082015-06-23T10:35:46.488-04:002015-06-23T10:35:46.488-04:00Boyd's book is free to download in pdf or brow...Boyd's book is free to <a href="https://scholar.google.com/scholar?q=boyd+chebyshev+and+fourier+spectral+methods" rel="nofollow">download in pdf</a> or <a href="https://books.google.com/books?id=lEWnQWyzLQYC&lpg=PA1&ots=WSdLuFo4uu&dq=boyd%20chebyshev%20and%20fourier%20spectral%20methods&lr&pg=PA1#v=onepage&q&f=false" rel="nofollow">browse on google books</a>.Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-60160810302659077212015-06-23T10:24:10.298-04:002015-06-23T10:24:10.298-04:00See the sections in Chapter 6 of Boyd's book w...See the sections in Chapter 6 of Boyd's book with the heading "Special Problems in Higher Dimensions:" for some pointers that might be helpful. Good luck with your calculation; would love to see your results when you get them. Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-84703092080092852722015-06-23T09:52:46.719-04:002015-06-23T09:52:46.719-04:00Hi Amit,
It's hard to tell what might be goi...Hi Amit, <br />It's hard to tell what might be going wrong with your code without more details. What grid point distribution are you using? For your periodic tests you probably used evenly spaced points. Are you using a Chebyshev roots grid for your non-periodic tests? Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-23781891628291629442015-06-02T22:26:36.083-04:002015-06-02T22:26:36.083-04:00Great blog! I am trying to compute laplacian of a...Great blog! I am trying to compute laplacian of a 2D function (non-periodic) using FFT technique. I first tested periodic functions and things work just fine. I am trying to extent your example to 2D function but cannot seem to figure our the recurrence for the derivative. My trial function is exp(x+y) sampled on grid with Nx and Ny points along x and y axis. The analytical result is 2.0*exp(x+y). I can send you my code if you want to see - but I have no luck so far. Any clues? Please help. Thx AmitAnonymoushttps://www.blogger.com/profile/17215657575191151707noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-56681724777334725252012-12-06T20:12:24.432-05:002012-12-06T20:12:24.432-05:00Great Blog... thanks for helping out the lot of us...Great Blog... thanks for helping out the lot of us wandering about details of implementation (and guts) of these methods... Keep posting, please!!! Unterschriberhttps://www.blogger.com/profile/07241636535437987992noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-77324220289133356772011-04-21T08:07:01.735-04:002011-04-21T08:07:01.735-04:00Thanks Anonymous; I filed a TOS complaint with Blo...Thanks Anonymous; I filed a TOS complaint with Blogger against that farm (it's funny, they even kept the links to my previous posts intact).Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-53772858779349264932011-04-21T05:31:00.331-04:002011-04-21T05:31:00.331-04:00You might be interested in the fact that your post...You might be interested in the fact that your post got copied over here. Looks like a farm anyway.<br /><br />http://onlinegeometryproblems.blogspot.com/2010/08/fftw-discrete-cosine-transform.html<br /><br />Oh and by the way: Very informative post - doing something similar myself atm.<br /><br />cheersAnonymousnoreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-68332803038940596662010-01-27T05:48:09.918-05:002010-01-27T05:48:09.918-05:00Hi jstults,
I found the following sentence in &qu...Hi jstults,<br /><br />I found the following sentence in "Spectral Methods" by C.Canuto in chapter 2.4.1:<br />"Virtually all of the classical literature on Chebyshev spectral methods uses this reversed order. Therefore, in the special case of the Chebyshev quadrature points we shall adhere to the ordering convention that is widely used int the literature (and implemented in the available software)."<br /><br />Greetings<br /> JensJenshttp://www.nourl.orgnoreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-38768000671271430792010-01-19T10:23:53.577-05:002010-01-19T10:23:53.577-05:00Thank you again Jens, it's nice to get all the...Thank you again Jens, it's nice to get all these little gotchas up for people to learn from. These down in the code details don't tend to make it into text books (but they tend to be a big stumbling block for learners).Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-1035662981222629162010-01-19T10:18:49.035-05:002010-01-19T10:18:49.035-05:00Hi,
another comment ;)
I added the fourier-fourie...Hi,<br /><br />another comment ;)<br />I added the fourier-fourier-chebyshev to p3dfft and finally found the real reason for your beta=-beta problem.<br />It is actually quite simple, but important to know:<br /><br />I used REDFT00 for which FFTW expects to chebyshev-nodes to be calculated with:<br />Fortran: coordZ(k) = cos(pi*(k-1)/(n3-1)) <br />C++ : coordZ(k) = cos(pi*k/double(n3))<br />=> coordZ needs to run from +1 to -1 (not -1 to +1)<br /><br />I suppose FFTW expects the coordinates to run from +1 to -1 for REDFC01, too.<br /><br />Greetings<br />JensJenshttp://www.nourl.orgnoreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-50468309880395598802010-01-15T08:41:57.824-05:002010-01-15T08:41:57.824-05:00p3dfft does the decomposition+3dfft of a 3d-grid o...p3dfft does the decomposition+3dfft of a 3d-grid over thousands of CPUs ...<br />internal it does the following:<br />state: all cpus have complete lines of real-data (phys.space) in x-direction<br />1) 1d-fft of all x-lines on all cpus (r2c)<br />2) rotate data that all cpus have complete lines of data in y-direction<br />3)1d-fft of all y-lines on all cpus (c2c)<br />4)rotate data that all cpus have complete lines of data in z-direction<br />5)1d-fft of all z-lines on all cpus (c2c)<br />=> 3d-fft with distributed memory (tested in my case with 16k cpus on BlueGene/P)<br /><br />I found a GPL-code for 3d-channelflow on www.channelflow.org.<br />It combines fourier+chebyshev+fourier and does an differentiation(REDCT00) using this transformed data.<br />The following functions describe all this:<br />FlowField.h<br /> void FlowField::makeSpectral()<br /> void FlowField::makePhysical()<br />diffops.h<br /> FlowField diff(..)<br /> FlowField xdiff(..)<br /> FlowField ydiff(..)<br /> FlowField zdiff(..)<br /><br />I would like to extend p3dfft with that feature.Jenshttp://www.nourl.orgnoreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-46574282613290271442010-01-15T07:51:45.541-05:002010-01-15T07:51:45.541-05:00Oh, I see, p3dfft is an MPI/Fortran wrapper around...Oh, I see, <a href="http://www.sdsc.edu/us/resources/p3dfft/index.php" rel="nofollow">p3dfft</a> is an MPI/Fortran wrapper around other FFT libraries like FFTW, nice.Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-19775360238937380992010-01-15T07:42:23.476-05:002010-01-15T07:42:23.476-05:00Jens: using the lib "p3dfft"
Why not fft...Jens: <i>using the lib "p3dfft"</i><br />Why not fftw? Is there something you especially like about 'pd3fft' (I'm always curious about exploring options)?<br /><br />I have been meaning to try a similar project for a flat plate boundary layer study (so it would be DCTs in stream-wise and wall-normal, and FFT's span-wise).<br /><br />I haven't actually done it yet though, so I couldn't give you any practical implementation suggestions or pitfalls. If I recall correctly (it's been a while) Boyd describes a mixed approach like this applied to a global circulation model <a href="http://books.google.com/books?id=lEWnQWyzLQYC&lpg=PP1&dq=boyd%20chebyshev%20and%20fourier%20spectral%20methods&pg=PP1#v=onepage&q=&f=false" rel="nofollow">in his book</a> (but doesn't go into much detail). If you find any good links please share, sorry I couldn't be of more help, and good luck!Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-47844913808076400252010-01-15T04:38:47.545-05:002010-01-15T04:38:47.545-05:00Hi jstults,
for computational fluid dynamics (CFD...Hi jstults,<br /><br />for computational fluid dynamics (CFD) I am calculating a lot on 3D-grids with periodic boundaries in all directions using the lib "p3dfft".<br />Calculating derivatives for this is easy, because you can do 3d-fft. Multiplying the 3d-result in fourier-space with a certain fixed 3d-array and transforming it back to physical space returns the derivative in x,y or z.<br /><br />Now my question:<br />I would like to change the boundaries that I only have periodicity in x and y but not in z (calculating channel flows). Therfor I want to use chebyshev-transform in z.<br />Is it possible to recive derivatives by forward-transform a 3d-array(physical space) to a 3d-array(fourier space) using r2r_cheby(in x) than r2c_fourier(in y) than c2c_fourier(in z). I wonder if it is "legal" to mix 2xDFT and 1xDCT and how to modify the resulting complex-data to get derivatives after backward-transform.Jenshttp://www.nourl.orgnoreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-84873423035496967692010-01-14T14:56:54.291-05:002010-01-14T14:56:54.291-05:00Good catch Jens, thank you.Good catch Jens, thank you.Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-13845725918308248142010-01-14T10:17:43.725-05:002010-01-14T10:17:43.725-05:00the last post also explains the "beta=-beta&q...the last post also explains the "beta=-beta"-problem you have:<br />because y[0..n-1]=-cos(pi/(2n) *(2j+1) this minus-sign should go into x of cos(n*arccos(x)) and therefor into Xj of the DCTs.<br /><br />Greetings<br />JensJenshttp://www.nourl.orgnoreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-75211750655425869642010-01-14T08:43:16.194-05:002010-01-14T08:43:16.194-05:00Hi Jens,
Glad it's helpful; you also may want...Hi Jens,<br /><br />Glad it's helpful; you also may want to look at <a href="http://j-stults.blogspot.com/2009/07/derivative-of-unequally-spaced-points.html" rel="nofollow">Calculating Derivatives of Unequally Spaced Points</a>, you don't <i>have</i> to use a Chebyshev-roots grid (but it is optimal).Joshua Stultshttps://www.blogger.com/profile/03506970399027046387noreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-89609762958840744262010-01-14T07:44:42.239-05:002010-01-14T07:44:42.239-05:00sorry... arb_int_opt_nodes(..) is correct and retu...sorry... arb_int_opt_nodes(..) is correct and returns<br />y[0..n-1] = cos(pi-pi/(2n)) .. cos(pi/(2n))<br /><br />----<br />Because REDFT10 is shifted by j+0.5 (j[0..n-1]) the chebyshev-nodes of the input data are set at<br />y[j=0..n-1] = +cos(pi -pi/n*(j+0.5))<br /> = -cos(pi/n*(j+0.5))<br /> = -cos(pi/(2n) *(2j+1))<br />You are using<br />y[n-i-1] = cos((2i+1)*pi/(2n)), i[0..n-1]<br />which is the same result.<br /><br />Greetings<br />JensJenshttp://www.gmx.denoreply@blogger.comtag:blogger.com,1999:blog-5822805028291837738.post-37057126226466470822010-01-14T05:36:49.155-05:002010-01-14T05:36:49.155-05:00Hi jstults,
thanks for your blog ... it really he...Hi jstults,<br /><br />thanks for your blog ... it really helps a lot.<br />I am still working on figuring out all details of chebyshev ...<br /><br />in arb_int_opt_nodes(..) you define the chebyshev-nodes at<br /> y[0..n-1] = cos(pi) .. cos(pi/(2n))<br />shouldn't that be at<br />y[0..n-1] = cos(pi-pi/(2n)) .. cos(pi/(2n)) for<br />REDFT10 ?<br /><br />Greetings<br />JensJenshttp://www.nourl.orgnoreply@blogger.com