|
[Sponsors] |
December 26, 2015, 10:20 |
Any experts on FFTW?
|
#1 |
New Member
Join Date: Dec 2015
Posts: 5
Rep Power: 11 |
Hi,
I recently got a Fortran77 DNS code and I would like to replace all the fft-replated routines there with FFTW3 (currently the code is equipped with FFTPACK. I notice that the code will blow up in the current setting; that's why I'm trying to use FFTW3.) By the way, is it known that FFTPACK is less stable than FFTW3? I have no idea about this. It's just my supervisor asked me to try that..... Now I'm playing with FFTW3. I am stuck at the call of dfftw_plan_many_dft_r2c_, which is a routine for transforming many real-data arrays of same size to complex. In Fortran, the layout of data structure is column-major. The following is what I have tried. I would like to column-transform or row transform a 2D real data array. For real data a(Nx,Nz), I tried first to transform the columns, that is Nz 1D array of size Nx. The plan is designed as follows call dfftw_plan_many_dft_r2c_(plan_forward,1,Nx,Nz, a,Nx,1,Nx, b,Nxc,1,Nxc,FFTW_ESTIMATE) where b is complex output data of size (Nxc,Nz) and Nxc=Nx/2+1. The results can be understood. But when I tried to transform the rows, that is Nx 1D array of size Nz. I design the plan as call dfftw_plan_many_dft_r2c_(plan_forward2,1,Nz,Nx, a,Nz,Nz,1, d,Nzc,Nzc,1,FFTW_ESTIMATE) where d is complex output data of size (Nx,Nzc) and Nzc=Nz/2+1. The results are wrong. I can't understand why that's the case. Is there any problem of my plan designing? Thanks for your help. jinhua |
|
December 26, 2015, 17:31 |
|
#2 |
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,290
Rep Power: 34 |
long time ago i used fftw_plan_dft_1d
if you want i could share that peace of code if you give me your email. r2c was also not very different to use if i rememebr correctly. It was 2010 to 2011 so i dont remember well. |
|
December 28, 2015, 11:07 |
|
#3 |
Member
Kaya Onur Dag
Join Date: Apr 2013
Posts: 94
Rep Power: 13 |
I have been there some time ago. Below some routines that I used to be able to understand how this library works specifically for taking derivatives (so you can also see the wavenumber order). Let me know if you can't figure out my messy coding, or still have the issue.
fftw 1d derivative Code:
program try_1d_diff use, intrinsic :: iso_c_binding implicit none include 'fftw3.f03' TYPE(C_PTR) :: pf,pb real*8 :: z REAL*8 , PARAMETER :: pi=4*ATAN(1.0D0) INTEGER, PARAMETER :: Nz=1*1e6 REAL*8, PARAMETER :: dz=2.0D0*pi/Nz COMPLEX(C_DOUBLE_COMPLEX), dimension(Nz) :: w,wh, dwdzh, dwdz INTEGER(C_INT):: i,kz(Nz) COMPLEX(C_DOUBLE_COMPLEX), parameter :: i_ = (0.0D0,1.0D0) kz = (/(I, I = 0, Nz/2-1), (I, I = -Nz/2,-1,1)/) !print*, 'kz:',kz; !print*, '' pf = fftw_plan_dft_1d(Nz,w,wh,-1,FFTW_ESTIMATE) pb = fftw_plan_dft_1d(Nz,dwdzh,dwdz,1,FFTW_ESTIMATE) do i = 1,Nz z = (i-1.0D0)*dz w(i) = sin(4*z) + cos(4*z)*i_ end do CALL FFTW_EXECUTE_DFT(pf,w,wh) print*,wh(5) !do i=1,Nz ! dwdzh(i)=i_*kz(i)*wh(i) !end do ! get back to the real space CALL FFTW_EXECUTE_DFT(pb,dwdzh,dwdz) !print*, 'wh', wh; ! print*, '' !print*, 'dwdzh', dwdzh; ! print*, '' !print*, 'dwdz', dwdz/8; !print*, ' ' CALL FFTW_DESTROY_PLAN(pf) CALL FFTW_DESTROY_PLAN(pb) end program try_1d_diff !!$ ifort -fr try_fft.f90 -I ~/fortran_lib/fftw/include -L fortran_lib/fftw/lib -lfftw3 -o pgm Code:
program try_2d_complex_diff use, intrinsic :: iso_c_binding implicit none include 'fftw3.f03' TYPE(C_PTR) :: pfft, pifft INTEGER(C_INT), PARAMETER :: Nx = 60 , Ny = 80 REAL(C_DOUBLE) :: dx,dy,pi = 4*atan(1.0D0),errx(Ny,Nx),erry(Ny,Nx) COMPLEX(C_DOUBLE_COMPLEX), dimension(Ny,Nx) :: u,dx_u,dy_u, u_h, dx_u_h, dy_u_h INTEGER(C_INT):: i,j,ky(Ny,Nx),kx(Ny,Nx) COMPLEX(C_DOUBLE_COMPLEX), parameter :: i_ = (0.0D0,1.0D0) dx = 2*pi/Nx dy = 2*pi/Ny do i=0,Nx-1 do j=0,Ny-1 u(j+1,i+1)=(exp(sin(1.0D0*i*dx))*exp(cos(1.0D0*j*dy))) end do end do print *, u print*, ' ----------- ' kx(1,:) = (/(I, I = 0, Nx/2), (I, I = -Nx/2+1,-1,1)/) ky(:,1) = (/(I, I = 0, Ny/2), (I, I = -Ny/2+1,-1,1)/) do i = 1,Nx kx(:,i)=kx(1,i) end do do j = 1,Ny ky(j,:)=ky(j,1) end do !!$ write(*,*) 'kx:' !!$ do j = 1,Ny !!$ write(*,*) kx(j,:) !!$ end do !!$ write(*,*) 'ky:' !!$ do j = 1,Ny !!$ write(*,*) ky(j,:) !!$ end do pfft = fftw_plan_dft_2d(Nx,Ny,u,u_h,-1,FFTW_ESTIMATE) !plan fft pifft = fftw_plan_dft_2d(Nx,Ny,dx_u_h,dx_u,1,FFTW_ESTIMATE) !plan ifft CALL FFTW_EXECUTE_DFT(pfft,u,u_h) dx_u_h=i_*kx*u_h dy_u_h=i_*ky*u_h CALL FFTW_EXECUTE_DFT(pifft,dx_u_h,dx_u) CALL FFTW_EXECUTE_DFT(pifft,dy_u_h,dy_u) print*, '------------' ! check dif y do i=1,Nx do j=1,Ny ! print*, dx_u(j,i)/Nx/Ny erry(j,i) =abs( dy_u(j,i)/Nx/Ny - (-exp(cos((j-1.0D0)*dy))*exp(sin((i-1.0D0)*dx))& &*sin((j-1.0D0)*dy))) errx(j,i) =abs( dx_u(j,i)/Nx/Ny - (exp(cos((j-1)*dy))& &*exp(sin((i-1)*dx))*cos((i-1)*dx))) end do end do print*,'errx:', maxval(errx) print*,'erry:', maxval(erry) CALL FFTW_DESTROY_PLAN(pfft) CALL FFTW_DESTROY_PLAN(pifft) end program try_2d_complex_diff !!$ ifort -fr try_fft.f90 -I ~/fortran_lib/fftw/include -L fortran_lib/fftw/lib -lfftw3 -o pgm Code:
program try_2d_diff use, intrinsic :: iso_c_binding implicit none include 'fftw3.f03' TYPE(C_PTR) :: pfft, pifft INTEGER(C_INT), PARAMETER :: Nx = 20 , Ny = 20 REAL(C_DOUBLE), dimension(Ny,Nx) :: u,dx_u,dy_u REAL(C_DOUBLE) :: dx,dy,pi = 4*atan(1.0D0), errx(Ny,Nx),erry(Ny,Nx)& &,t,t1,t2 COMPLEX(C_DOUBLE_COMPLEX), dimension(Ny/2+1,Nx) :: u_h, dx_u_h, dy_u_h INTEGER(C_INT):: i,j,ky(Ny/2+1),kx(Nx) COMPLEX(C_DOUBLE_COMPLEX), parameter :: i_ = (0.0D0,1.0D0) dx = 2*pi/Nx dy = 2*pi/Ny ! print*, u ! print*, ' ----------- ' kx = (/(I, I = 0, Nx/2), (I, I = -Nx/2+1,-1,1)/) ky = (/(I, I = 0, Ny/2)/) ! print*, 'kx:' , kx,'ky:',ky pfft = fftw_plan_dft_r2c_2d(Nx,Ny,u,u_h,FFTW_MEASURE) !plan fft pifft = fftw_plan_dft_c2r_2d(Nx,Ny,dx_u_h,dx_u,FFTW_MEASURE) !plan !ifft do i=0,Nx-1 do j=0,Ny-1 u(j+1,i+1)=exp(sin(1.0D0*i*dx))*exp(cos(1.0D0*j*dy)) end do end do call cpu_time(t1) CALL FFTW_EXECUTE_DFT_r2c(pfft,u,u_h) call cpu_time(t2) do i = 1,Nx do j = 1,Ny/2+1 dx_u_h(j,i)=i_*kx(i)*u_h(j,i) dy_u_h(j,i)=i_*ky(j)*u_h(j,i) end do end do t=t2-t1 call cpu_time(t1) CALL FFTW_EXECUTE_DFT_c2r(pifft,dx_u_h,dx_u) call cpu_time(t2) t=t+t2-t1 CALL FFTW_EXECUTE_DFT_c2r(pifft,dy_u_h,dy_u) print*, 'time:', t ! check dif y do i=1,Nx do j=1,Ny ! print*, dx_u(j,i)/Nx/Ny erry(j,i) =abs( dy_u(j,i)/Nx/Ny - (-exp(cos((j-1.0D0)*dy))*exp(sin((i-1.0D0)*dx))& &*sin((j-1.0D0)*dy))) errx(j,i) =abs( dx_u(j,i)/Nx/Ny - (exp(cos((j-1)*dy))& &*exp(sin((i-1)*dx))*cos((i-1)*dx))) end do end do print*,'errx:', maxval(errx) print*,'erry:', maxval(erry) CALL FFTW_DESTROY_PLAN(pfft) CALL FFTW_DESTROY_PLAN(pifft) end program try_2d_diff !!$ ifort -fr try_fft.f90 -I ~/fortran_lib/fftw/include -L fortran_lib/fftw/lib -lfftw3 -o pgm Code:
program try_2d_diff use, intrinsic :: iso_c_binding implicit none include 'fftw3.f03' TYPE(C_PTR) :: pfft, pifft INTEGER(C_INT), PARAMETER :: Nx = 1024 , Ny = 1024 REAL(C_DOUBLE), dimension(Ny,Nx) :: u,dx_u,dy_u REAL(C_DOUBLE) :: dx,dy,pi = 4*atan(1.0D0), errx(Ny,Nx),erry(Ny,Nx)& &,t,t1,t2 COMPLEX(C_DOUBLE_COMPLEX), dimension(Ny/2+1,Nx) :: u_h, dx_u_h, dy_u_h INTEGER(C_INT):: i,j,ky(Ny/2+1),kx(Nx),rank,howmany,idist& &,odist,istride,ostride,inembed(2),onembed(2),n(2),n_h(2) COMPLEX(C_DOUBLE_COMPLEX), parameter :: i_ = (0.0D0,1.0D0) rank = 2 !Two dimensional howmany = 1 !2D, one-shot fft n = (/Nx,Ny/) !input- slow one first comes n_h = (/Nx,Ny/2+1/) !dimension in fourier space idist = 0 !unused because how many is one odist = 0 !unused because how many is one istride = 1 !array is contagious in the memory ostride = 1 !array is contagious in the memory inembed = n !entering dimesions onembed = n_h !exiting dimesions pfft = fftw_plan_many_dft_r2c(rank,n,howmany,u,inembed,istride& &,idist,u_h,onembed,ostride,odist,FFTW_MEASURE); ! forward transform pifft = fftw_plan_many_dft_c2r(rank,n,howmany,dx_u_h,onembed,ostride& &,odist,dx_u,inembed,istride,idist,FFTW_MEASURE); ! forward transform !calculate wavenumbers kx = (/(I, I = 0, Nx/2), (I, I = -Nx/2+1,-1,1)/) ky = (/(I, I = 0, Ny/2)/) !pfft = fftw_plan_dft_r2c_2d(Nx,Ny,u,u_h,FFTW_MEASURE) !plan fft !pifft = fftw_plan_dft_c2r_2d(Nx,Ny,dx_u_h,dx_u,FFTW_MEASURE) !plan !generate the mesh and function values dx = 2*pi/Nx dy = 2*pi/Ny do i=0,Nx-1 do j=0,Ny-1 u(j+1,i+1)=exp(sin(1.0D0*i*dx))*exp(cos(1.0D0*j*dy)) end do end do call cpu_time(t1) CALL FFTW_EXECUTE_DFT_r2c(pfft,u,u_h) call cpu_time(t2) do i = 1,Nx do j = 1,Ny/2+1 dx_u_h(j,i)=i_*kx(i)*u_h(j,i) dy_u_h(j,i)=i_*ky(j)*u_h(j,i) end do end do t=t2-t1 call cpu_time(t1) CALL FFTW_EXECUTE_DFT_c2r(pifft,dx_u_h,dx_u) call cpu_time(t2) t=t+t2-t1 CALL FFTW_EXECUTE_DFT_c2r(pifft,dy_u_h,dy_u) print*, 'time:', t ! check dif y do i=1,Nx do j=1,Ny ! print*, dx_u(j,i)/Nx/Ny erry(j,i) =abs( dy_u(j,i)/Nx/Ny - (-exp(cos((j-1.0D0)*dy))*exp(sin((i-1.0D0)*dx))& &*sin((j-1.0D0)*dy))) errx(j,i) =abs( dx_u(j,i)/Nx/Ny - (exp(cos((j-1)*dy))& &*exp(sin((i-1)*dx))*cos((i-1)*dx))) end do end do print*,'errx:', maxval(errx) print*,'erry:', maxval(erry) CALL FFTW_DESTROY_PLAN(pfft) CALL FFTW_DESTROY_PLAN(pifft) end program try_2d_diff !!$ ifort -fr try_fft.f90 -I ~/fortran_lib/fftw/include -L fortran_lib/fftw/lib -lfftw3 -o pgm Code:
program try_2d_diff_many_malloc use, intrinsic :: iso_c_binding implicit none include 'fftw3.f03' TYPE(C_PTR) :: pfft, pifft INTEGER(C_INT), PARAMETER :: Nx = 8 , Ny = 12 real(C_DOUBLE), pointer :: u(:,:)=>NULL(),dx_u(:,:)=>NULL(),dy_u(:,:)=>NULL() integer(C_SIZE_T) :: sz_real = Nx*Ny, sz_complex = Nx*(Ny/2+1) type(C_PTR) :: p_u,p_dx_u,p_dy_u,p_u_h,p_dx_u_h,p_dy_u_h ! REAL(C_DOUBLE), dimension(Ny,Nx) :: u,dx_u,dy_u REAL(C_DOUBLE) :: dx,dy,pi = 4*atan(1.0D0), errx(Ny,Nx),erry(Ny,Nx),t,t1,t2 complex(C_DOUBLE_COMPLEX), pointer :: u_h(:,:)=>null(),dx_u_h(:,:)=>null(),dy_u_h(:,:)=>null() ! COMPLEX(C_DOUBLE_COMPLEX), dimension(Ny/2+1,Nx) :: u_h, dx_u_h, dy_u_h INTEGER(C_INT):: i,j,ky(Ny/2+1),kx(Nx),rank,howmany,idist,odist,istride,ostride,inembed(2),onembed(2),n(2),n_h(2) COMPLEX(C_DOUBLE_COMPLEX), parameter :: i_ = (0.0D0,1.0D0) p_u=fftw_alloc_real(sz_real) p_dx_u = fftw_alloc_real(sz_real) p_dy_u = fftw_alloc_real(sz_real) p_u_h = fftw_alloc_complex(sz_complex) p_dx_u_h = fftw_alloc_complex(sz_complex) p_dy_u_h = fftw_alloc_complex(sz_complex) call c_f_pointer(p_u, u, [2:13,Nx ]) call c_f_pointer(p_dx_u, dx_u, [Ny,Nx ]) call c_f_pointer(p_dy_u, dy_u, [Ny,Nx ]) call c_f_pointer(p_u_h, u_h, [Ny/2+1,Nx]) call c_f_pointer(p_dx_u_h,dx_u_h,[Ny/2+1,Nx]) call c_f_pointer(p_dy_u_h,dy_u_h,[Ny/2+1,Nx]) print*, 'here' print*,'Ny:', LBOUND(u,1),UBOUND(u,1) print*,'Nx:', LBOUND(u,2),UBOUND(u,2) rank = 2 !Two dimensional howmany = 1 !2D, one-shot fft n = (/Nx,Ny/) !input- slow one first comes n_h = (/Nx,Ny/2+1/) !dimension in fourier space idist = 0 !unused because how many is one odist = 0 !unused because how many is one istride = 1 !array is contagious in the memory ostride = 1 !array is contagious in the memory inembed = n !entering dimesions onembed = n_h !exiting dimesions pfft = fftw_plan_many_dft_r2c(rank,n,howmany,u,inembed,istride,idist,u_h,onembed,ostride,odist,FFTW_MEASURE); ! forward ! transform pifft = fftw_plan_many_dft_c2r(rank,n,howmany,dx_u_h,onembed,ostride& &,odist,dx_u,inembed,istride,idist,FFTW_MEASURE); ! forward transform !!$ !calculate wavenumbers kx = (/(I, I = 0, Nx/2), (I, I = -Nx/2+1,-1,1)/) ky = (/(I, I = 0, Ny/2)/) ! pfft = fftw_plan_dft_r2c_2d(Nx,Ny,u,u_h,FFTW_MEASURE) !plan fft ! pifft = fftw_plan_dft_c2r_2d(Nx,Ny,dx_u_h,dx_u,FFTW_MEASURE) !plan !generate the mesh and function values dx = 2*pi/Nx dy = 2*pi/Ny do i=0,Nx-1 do j=0,Ny-1 u(j+1,i+1)=(sin(4.0D0*i*dx)) end do end do call cpu_time(t1) CALL FFTW_EXECUTE_DFT_r2c(pfft,u,u_h) call cpu_time(t2) do i = 1,Nx do j = 1,Ny/2+1 dx_u_h(j,i)=i_*kx(i)*u_h(j,i) dy_u_h(j,i)=i_*ky(j)*u_h(j,i) end do end do t=t2-t1 call cpu_time(t1) CALL FFTW_EXECUTE_DFT_c2r(pifft,dx_u_h,dx_u) call cpu_time(t2) t=t+t2-t1 CALL FFTW_EXECUTE_DFT_c2r(pifft,dy_u_h,dy_u) print*, 'time:', t ! check dif y do i=1,Nx do j=1,Ny ! print*, dx_u(j,i)/Nx/Ny erry(j,i) =abs( dy_u(j,i)/Nx/Ny - 0 ) errx(j,i) =abs( dx_u(j,i)/Nx/Ny - 4*(cos(4.0D0*i*dx))) end do end do print*,'errx:', maxval(errx) print*,'erry:', maxval(erry) CALL FFTW_DESTROY_PLAN(pfft) CALL FFTW_DESTROY_PLAN(pifft) call fftw_free(p_u) call fftw_free(p_dx_u) call fftw_free(p_dy_u) call fftw_free(p_u_h) call fftw_free(p_dx_u_h) call fftw_free(p_dy_u_h) end program try_2d_diff_many_malloc !!$ifort -fr try_fft.f90 -I ~/fortran_lib/fftw/include -L fortran_lib/fftw/lib -lfftw3 -o pgm -I suggest you to use 'many' and 'malloc' way of doing this. happy new year! |
|
May 22, 2022, 20:44 |
Similar Problem
|
#4 | |
New Member
Negar
Join Date: May 2022
Posts: 2
Rep Power: 0 |
Quote:
I have a similar issue and could not solve it even after searching for so long! I was wondering if your problem is solved. Would you mind taking a look at my thread? Partial Derivative Using Fourier Transform (FFTW) in 2D I am using FFTW3 in C++ for derivation. Many thanks, Negar |
||
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Experts in Simulating BUbble Column Reactors | aliyah | Main CFD Forum | 7 | June 19, 2012 03:38 |
Post-processing U field with fftw | Pascal_doran | OpenFOAM Post-Processing | 1 | November 23, 2010 04:05 |
FFTW for Fortran | Benjamin Cassart | Main CFD Forum | 3 | August 1, 2007 11:31 |
help me experts | chelesa | CFX | 0 | February 15, 2004 04:44 |
invitation for CFD experts | ROOH-UL-AMIN KHURRAM | Main CFD Forum | 2 | March 17, 2000 01:02 |