CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Any experts on FFTW?

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   December 26, 2015, 10:20
Default Any experts on FFTW?
  #1
New Member
 
Join Date: Dec 2015
Posts: 5
Rep Power: 8
jinhua2015 is on a distinguished road
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
jinhua2015 is offline   Reply With Quote

Old   December 26, 2015, 17:31
Default
  #2
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,053
Rep Power: 28
arjun will become famous soon enougharjun will become famous soon enough
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.
arjun is offline   Reply With Quote

Old   December 28, 2015, 11:07
Default
  #3
Member
 
Kaya Onur Dag
Join Date: Apr 2013
Posts: 94
Rep Power: 10
kaya is on a distinguished road
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
fftw 2d complex2complex derivative
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
fftw 2d complex2real
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
the same with using 'many' which allows you to take 2dfft on 3d data without using a loop
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
with c malloc allocations - as suggested in the fftw manual
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
-they might not work if you directly try to run them but you should be able to get what you actually need from the coding.
-I suggest you to use 'many' and 'malloc' way of doing this.

happy new year!
kaya is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 12:45.