Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Topic: using FFTW in physics
Replies: 0

 roman Posts: 4 Registered: 7/16/09
using FFTW in physics
Posted: Jul 19, 2013 7:30 PM

Dear all.

I am writing a simulation code that uses FFTW. I have not used FFTW before and thus I want to make sure that I understand some things correctly. Let's say we have a 3D real-valued field fr that has nx,ny,nz components. I work with Fortran, so the declaration is fr(nx,ny,nz). FFTW gives me a complex-valued Fourier image of this field fk that has nx/2+1,ny,nz components. This is how I do it:

double precision :: fr(nx,ny,nz)
double complex :: fk(nx/2+1,ny,nz)

call dfftw_plan_dft_r2c_3d(plan, nx, ny, nz, fr, fk, FFTW_ESTIMATE)
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)

The inverse transform is done this way:

call dfftw_plan_dft_c2r_3d(plan, nx, ny, nz, fk, fr, FFTW_ESTIMATE)
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)
fr = fr/dble(nx*ny*nz)

Up to this, I am confident that it is correct. Now, how are the wavevectors assigned to the individual values in fk? I know that fk(1,1,1) corresponds to the dc component and all indices of fk represent positive wavevectors. Therefore, if the size of my simulated box in the real space is Lx=nx*ax, Ly=ny*ay, Lz=nz*az (ax,ay,az are the lattice parameters along the three axes), the indices of fk should correspond to these wavevector components:

1st index of fk: 2*pi/Lx*(0, ..., nx/2) = (0, ..., pi/ax)
2nd index of fk: 2*pi/Ly*(0, ..., ny-1) = (0, ..., 2*pi/ay - 2*pi/Ly)
3rd index of fk: 2*pi/Lz*(0, ..., nz-1) = (0, ..., 2*pi/az - 2*pi/Lz)

Do I understand the assignment of the wavevectors above to individual indices of fk correctly? I will eventually have to differentiate fk w.r.t. kx,ky,kz and so this is very important.