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 realvalued field fr that has nx,ny,nz components. I work with Fortran, so the declaration is fr(nx,ny,nz). FFTW gives me a complexvalued 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, ..., ny1) = (0, ..., 2*pi/ay  2*pi/Ly) 3rd index of fk: 2*pi/Lz*(0, ..., nz1) = (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.
Thanks a lot for your comments,
Roman

