Reputation: 1
What am I missing from FFTW's advanced interface to make it the same as the basic interface for a test case?
I am using a legacy CFD code that uses a pseudospectral solver, and ultimately I want to update the code with better FFT algorithms, so I'm currently learning the interface of FFTW3 in Fortran. I have a 3D velocity array from the DNS code in various forms to test the FFTs and make sure I'm transforming things correctly.
Currently, I have a 3D array which is physical in the first dimension but spectral in the other two dimensions. My goal is to do a 2D transform on this array to get the 3D physical values, and I'm stuck trying to figure out what is wrong with my current code.
I have successfully done a 2D transform of one single plane of my 3D input array using fftw_plan_dft_c2r(), but when I tried to iterate over the planes, I got bad results (I assume because of something happening with the destroyed input vector, but when I tried to preserve the input, I got a segmentation fault).
I figured I might be better off using the advanced interface to do multiple DFTs on the planes, using fftw_plan_many_dft_c2r(). So the first test case is what I'm showing in the code provided, and that's doing a transform of a single 2D plane with both methods. For certain parameters, I believe these routines should be identical. The expected output is an mz x mx array of a constant value around 150. The basic interface routine indeed gives me this result, but the advanced interface gives me something that looks like
9576957.89326171 -6100449.91245205 4941.42749753309
-2029633.19975117 378.904251844637 -1219266.13730242
1476.66846120109 -868929.573548530 -20.9026499204410
-676593.055185564 450.197926313034 -552121.186219578
-425.926031798779 -467648.354118945 -220.256552425897
...
I'm really not sure what else I could even change. Out of curiosity, I played around with {i,o}stride, {i,o}dist, and {i,o}nembed, but nothing seemed to get me closer to the right answer.
Code below:
program fftw_test
use, intrinsic :: iso_c_binding
implicit none
include 'fftw3.f03'
! Grid sizes
integer, parameter :: nxh = 128, nz = 128
integer, parameter :: mx = 384, mz = 192
! FFTW Variables
type(C_PTR) :: plan2,plan3
complex(C_DOUBLE_COMPLEX), dimension(nz,nxh) :: u_phys
complex(C_DOUBLE_COMPLEX), dimension(nz,nxh) :: u1_2d,u3_2d
real(C_DOUBLE), dimension(mz,mx) :: u2_2d,u4_2d
integer(C_INT), dimension(2) :: msize_2d = (mx,mz)
integer(C_INT) :: rank, howmany, idist, odist, istride, ostride
! ------------------------------------------------------------------------------------ !
! ------------------------------------------------------------------------------------ !
! Create FFT plan (before initializing variable)
print *,'Create FFT Plan(s)'
rank = 2; howmany = 1; ! 1 2D c2r transform on an (mz) x (mx) array
! (msize_2d is flipped because of column-major arrays)
idist = 0; odist = 0; ! Unused since howmany = 1
istride = 1; ostride = 1; ! Array is contiguous in memory (I think?)
plan2 = fftw_plan_dft_c2r (rank, msize_2d,u1_2d,u2_2d,FFTW_ESTIMATE)
plan3 = fftw_plan_many_dft_c2r(rank, msize_2d,howmany, &
u3_2d,msize_2d,istride,idist, &
u4_2d,msize_2d,ostride,odist, FFTW_ESTIMATE)
u_phys = (0.0,0.0)
u_phys(1,1) = (150.569261744544,1.136868377216160E-013)
u1_2d = u_phys ! temp storage
u3_2d = u_phys ! temp storage
! Execute FFTs
print *,'Execute FFTs'
call fftw_execute_dft_c2r(plan2,u1_2d,u2_2d)
write(200,*) u2_2d
call fftw_execute_dft_c2r(plan3,u3_2d,u4_2d)
write(201,*) u4_2d
! Destroy FFT plans
call fftw_destroy_plan(plan2)
call fftw_destroy_plan(plan3)
end program fftw_test
EDIT: I updated the code to explicitly define the grid size variables and the sample spectral array to be transformed. However, for some reason the 2D transform no longer works for me despite not changing anything in the code from my previous verification.
Upvotes: 0
Views: 126
Answers (1)
Reputation: 1
For completeness, I'll share my solution to this problem as an answer. The following is a slightly altered code that indeed gives the expected output:
program fftw_test
use,intrinsic :: iso_c_binding
implicit none
include 'fftw3.f03'
! Grid sizes
integer,parameter :: nz = 128, nxh = 128
integer,parameter :: mz = 192, mx = 384
! FFTW Plans
type(C_PTR) :: plan1,plan2
! Complex arrays
complex(C_DOUBLE_COMPLEX), dimension(nz,nxh) :: u_phys ! input array
complex(C_DOUBLE_COMPLEX), dimension(mz,mx) :: u1_spec_2D
complex(C_DOUBLE_COMPLEX), dimension(mz,mx) :: u2_spec_2D
! Real arrays
real(C_DOUBLE), dimension(mz,mx) :: u1_real_2D ! 2D c2r transform target array
real(C_DOUBLE), dimension(mz,mx) :: u2_real_2D ! 2D c2r transform target array
! fft planl vars
integer(C_INT) :: rank, howmany, idist, odist, istride, ostride
integer(C_INT), dimension(2) :: msz = (/mz,mx/)
! ------------------------------------------------------------------------------------ !
! ------------------------------------------------------------------------------------ !
! Create FFT plan
rank = 2; howmany = 1
idist = 0; odist = 0
istride = 1; ostride = 1
plan1 = fftw_plan_dft_c2r(rank,msz,u1_spec_2D,u1_real_2D,FFTW_ESTIMATE)
plan2 = fftw_plan_many_dft_c2r(rank,msz,howmany, &
u2_spec_2D,msz,istride,idist, &
u2_real_2D,msz,ostride,odist,FFTW_ESTIMATE)
u_phys = (0.0,0.0)
u_phys(1,1) = (150.569261744544,1.136868377216160E-013)
u1_spec_2D(1:nz,1:nxh) = u_phys
u2_spec_2D(1:nz,1:nxh) = u_phys
! Execute FFT
call fftw_execute_dft_c2r(plan1,u1_spec_2D,u1_real_2D)
call fftw_execute_dft_c2r(plan2,u2_spec_2D,u2_real_2D)
write(100,*) u1_real_2D
write(200,*) u2_real_2D
! Destroy plans
call fftw_destroy_plan(plan1)
call fftw_destroy_plan(plan2)
end program fftw_test
Other than some name changes, the only difference is that I add extra space in the complex array to be transformed nz x nxh to mz x mx and I correctly make the size array an array instead of a complex number.
Upvotes: 0