FFTW based Two Dimensional Serial Poisson Solver

Program poisson2d

implicit none

integer ( kind = 4 ), parameter :: Nx = 128

integer ( kind = 4 ), parameter :: Ny = 128

integer ( kind = 4 ), parameter :: Nh = ( Nx / 2 ) + 1

real ( kind = 8 ), parameter :: pi=3.14159265358979323846d0

include "fftw3.f"

integer ( kind = 4 ) i,j

real ( kind = 8 ) Lx,Ly,dx,dy,kx,ky

real ( kind = 8 ) x(Nx),y(Ny),phi(Nx,Ny),rho(Nx,Ny),rho_dum(Nx,Ny)

complex ( kind = 8 ) rhok(Nh,Ny),rhok_dum(Nh,Ny),phik(Nh,Ny),phik_dum(Nh,Ny)

integer ( kind = 8 ) plan_forward,plan_backward

open(unit=10,file='poisson.dat',status='unknown')

Lx = 2.0*pi

Ly = 2.0*pi

dx = Lx/float(Nx)

dy = Ly/float(Ny)

do i = 1, Nx

x(i)=0.0d0+real(i-1)*dx

  do j = 1, Ny

  y(j)=0.0+real(j-1)*dy

  rho(i,j) = 2.0d0*sin(x(i))*cos(y(j))

  rho_dum(i,j) = rho(i,j)

  end do

end do

  call dfftw_plan_dft_r2c_2d_ (plan_forward, Nx, Ny, rho_dum, rhok, FFTW_ESTIMATE)

  call dfftw_execute_ (plan_forward)

do i = 1,Nx/2+1

  do j = 1,Ny/2

  kx = 2.0d0*pi*float(i-1)/Lx

  ky = 2.0d0*pi*float(j-1)/Ly

    if (i == 1 .and. j == 1) then

    phik(i,j) = (0.0d0,0.0d0)

    else

    phik(i,j) = rhok(i,j)/( kx*kx + ky*ky )

    endif

  phik_dum(i,j) = phik(i,j) 

  enddo

  do j = Ny/2+1,Ny

  kx = 2.0d0*pi*float(i-1)/Lx

  ky = 2.0d0*pi*float((j-1)-Ny)/Ly

  phik(i,j) = rhok(i,j)/( kx*kx + ky*ky )

  phik_dum(i,j) = phik(i,j) 

  enddo

enddo

  call dfftw_plan_dft_c2r_2d_ (plan_backward, Nx, Ny, phik_dum, phi, FFTW_ESTIMATE)

  call dfftw_execute_ (plan_backward)

  call dfftw_destroy_plan_ (plan_forward)

  call dfftw_destroy_plan_ (plan_backward)

do i = 1, Nx

  do j = 1, Ny

  phi(i,j) = phi(i,j)/(float(Nx)*float(Ny))

  write(10,*) x(i),y(j),rho(i,j),phi(i,j)

  end do

end do

end program poisson2d