Correct usage of zfft2d?

[abcduncan]

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.

Comer Duncan
Department of Physics and Astronomy
Bowling Green State University
Bowling Green, OH 43403
email: gcomerd@yahoo.com

!!!!!!!!!!!!!!!!! begin code fragment !!!!!!!!!!!!!!!
double precision,dimension(1:40000)::vpx,vpy,filxn2,filyn2,filxn1_0, &
filyn1_0,filxn2_0,filyn2_0,vpx_0,vpy_0
double precision,dimension(1:256,1:4096)::ux,uy,ffx,ffy,sx,sy,ss,aa,vvx,vvy,p,ux0,uy0
double complex,dimension(1:256,1:4096)::csx,csy,css
double complex,dimension(1:256,1:4096)::vxc,vyc,pc

!! snip snip
!

nx = 256
ny = 4096

j=1,ny
i=1,nx

sx(i,j)=dsin(2.d0*pi/dfloat(nx)*dfloat(i-1))

sy(i,j)=dsin(2.d0*pi/dfloat(ny)*dfloat(j-1))

ss(i,j)=sx(i,j)*sx(i,j)+sy(i,j)*sy(i,j)



aa(i,j)=1.d0+4.0d0*mu*(dt/(ro*h*h))*(dsin(pi*dfloat(i-1)/dfloat(nx))**2+ &
dsin(pi*dfloat(j-1)/dfloat(ny))**2)

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))

if(ss(i,j)==0.0d0) then

sx(i,j) = 0.0d0

sy(i,j) = 0.0d0

css(i,j) = 1.0d0

endif

enddo
enddo

!! snip snip snip

vxc=cmplx(vvx,0.d0)
vyc=cmplx(vvy,0.d0)

! call the complex to complex fft routines to get fft of vxc and vyc, the inputs to the solve
!

call zfft2d(vxc,nx,ny,-1)
call zfft2d(vyc,nx,ny,-1)

! now do arithmetic in k-space


pc=(sx*vxc+sy*vyc)/(css/2.d0)
vxc=(vxc-csx/2.d0*pc)/(aa/2.d0+0.5d0)
vyc=(vyc-csy/2.d0*pc)/(aa/2.d0+0.5d0)

! now do inverse fft to return new ux in vxc and new uy in vyc

call zfft2d(vxc,nx,ny,1)
call zfft2d(vyc,nx,ny,1)


! now assign real part of returned arrays to physical ux and uy
ux=real(vxc)
uy=real(vyc)
p =real(pc)


!!!!!!!!!!!!!! end code fragment !!!!!!!!!!!!!!!!!!!!!