Replacing numpy.fft routines with pyfftw, not working as expected

830 Views Asked by At

I have some working python code making use of the numpy.fft package, here is a snippet:

for i in range(steps):
    print i
    psixvec = Ux * psixvec
    psikvec = Uk * np.fft.fftn(psixvec)
    psixvec = np.fft.ifftn(psikvec)

return psixvec

I tried to rewrite this code to make use of the pyfftw package. What I came up with is the following code, which should work:

fft = fftw.builders.fftn(psix_align, auto_align_input = True, auto_contiguous = True,
                      overwrite_input = False, threads = 1, avoid_copy = False)

ifft = fftw.builders.ifftn(psik_align, auto_align_input = True, auto_contiguous = True,
                      overwrite_input = False, threads = 1, avoid_copy = False) 


for i in range(steps):
    psix_align[:] = Ux * psix_align
    psik_align[:] = Uk * fft()
    psix_align[:] = ifft()

return psix_align

The problem is, this code does not produce the same result as using the numpy.fft package. See the attached images.

numpy fft package pyfftw package

1

There are 1 best solutions below

0
On BEST ANSWER

Solved. For the initialization I was using

psix_align = fftw.n_byte_align(psi0, fftw.simd_alignment, dtype='complex64')
psik_align = fftw.n_byte_align(np.zeros_like(psi0), fftw.simd_alignment, dtype='complex64')

I needed to replace complex64 with complex128. Now I get the same result. That's probably due to the fact that the numbers involved are very small (see the 1e-11 on the z-axis).

Edit: Maybe someone can add pyfftw to the tags of the question?