I am currently working on a Navier-Stokes solver which assumes periodic boundary conditions in two dimensions. In the code, fragment shown below, I use the MKL zfft2d routine for complex to complex ffts. I am concerned that I may not be using it quite correctly and hope to get some advice on this from the community.
The basic issue is whether I need to use arrays of declared size 0:255 or 1:256. The NS solver in space uses 1:256 by 1:4096. I am getting some results which draw this into question. Also, I am not explicitly making the 1:256 by 1:4096 2D array periodic. I suspect that this can have bad consequences unless the fft routine actually takes care of assumed periodicity internally. I am also not calling the routine for any initialization prior to using it. From the documentation that I see, the initialization call sometimes needed by fft routines is not needed in this case. Correct?
To show you what I am doing I list below a code fragment and ask experienced users for their input. Thanks very much.
Department of Physics and Astronomy
Bowling Green State University
Bowling Green, OH 43403
!!!!!!!!!!!!!!!!! begin code fragment !!!!!!!!!!!!!!!
double precision,dimension(1:40000)::vpx,vpy,filxn2,filyn2,filxn1_0, &
!! snip snip
nx = 256
ny = 4096
csx(i,j) = cmplx(0.0d0,c*sx(i,j))
csy(i,j) = cmplx(0.0d0,c*sy(i,j))
css(i,j) = cmplx(0.0d0,c*ss(i,j))
sx(i,j) = 0.0d0
sy(i,j) = 0.0d0
css(i,j) = 1.0d0
!! snip snip snip
! call the complex to complex fft routines to get fft of vxc and vyc, the inputs to the solve
! now do arithmetic in k-space
! now do inverse fft to return new ux in vxc and new uy in vyc
! now assign real part of returned arrays to physical ux and uy
!!!!!!!!!!!!!! end code fragment !!!!!!!!!!!!!!!!!!!!!