FFTW based Three Dimensional Serial Poisson Solver
Program poisson3d
implicit none
integer ( kind = 4 ), parameter :: Nx = 64
integer ( kind = 4 ), parameter :: Ny = 64
integer ( kind = 4 ), parameter :: Nz = 64
integer ( kind = 4 ), parameter :: Nh = ( Nx / 2 ) + 1
real ( kind = 8 ), parameter :: pi=3.14159265358979323846d0
include "fftw3.f"
integer ( kind = 4 ) i,j,k
real ( kind = 8 ) Lx,Ly,Lz,dx,dy,dz,kx,ky,kz
real ( kind = 8 ) x(Nx),y(Ny),z(Nz),phi(Nx,Ny,Nz),rho(Nx,Ny,Nz),rho_dum(Nx,Ny,Nz)
complex ( kind = 8 ) rhok(Nh,Ny,Nz),rhok_dum(Nh,Ny,Nz),phik(Nh,Ny,Nz),phik_dum(Nh,Ny,Nz)
integer ( kind = 8 ) plan_forward,plan_backward
open(unit=10,file='poisson.dat',status='unknown')
Lx = 2.0*pi
Ly = 2.0*pi
Lz = 2.0*pi
dx = Lx/float(Nx)
dy = Ly/float(Ny)
dz = Ly/float(Nz)
do i = 1, Nx
x(i)=0.0d0+real(i-1)*dx
do j = 1, Ny
y(j)=0.0+real(j-1)*dy
do k = 1, Nz
z(k)=0.0+real(k-1)*dz
rho(i,j,k) = 3.0d0*sin(x(i))*cos(y(j))*sin(z(k))
rho_dum(i,j,k) = rho(i,j,k)
enddo
end do
end do
call dfftw_plan_dft_r2c_3d_ (plan_forward, Nx, Ny, Nz, rho_dum, rhok, FFTW_ESTIMATE)
call dfftw_execute_ (plan_forward)
do i = 1,Nh
do j = 1,Ny/2
do k = 1,Nz/2
kx = 2.0d0*pi*float(i-1)/Lx
ky = 2.0d0*pi*float(j-1)/Ly
kz = 2.0d0*pi*float(k-1)/Lz
if (i == 1 .and. j == 1 .and. k ==1) then
phik(i,j,k) = (0.0d0,0.0d0)
else
phik(i,j,k) = rhok(i,j,k)/( kx*kx + ky*ky + kz*kz)
endif
phik_dum(i,j,k) = phik(i,j,k)
enddo
do k = Nz/2+1,Nz
kx = 2.0d0*pi*float(i-1)/Lx
ky = 2.0d0*pi*float(j-1)/Ly
kz = 2.0d0*pi*float((k-1)-Nz)/Lz
phik(i,j,k) = rhok(i,j,k)/( kx*kx + ky*ky + kz*kz)
phik_dum(i,j,k) = phik(i,j,k)
enddo
enddo
do j = Ny/2+1,Ny
do k = 1,Nz/2
kx = 2.0d0*pi*float(i-1)/Lx
ky = 2.0d0*pi*float((j-1)-Ny)/Ly
kz = 2.0d0*pi*float(k-1)/Lz
if (i == 1 .and. j == 1 .and. k ==1) then
phik(i,j,k) = (0.0d0,0.0d0)
else
phik(i,j,k) = rhok(i,j,k)/( kx*kx + ky*ky + kz*kz)
endif
phik_dum(i,j,k) = phik(i,j,k)
enddo
do k = Nz/2+1,Nz
kx = 2.0d0*pi*float(i-1)/Lx
ky = 2.0d0*pi*float((j-1)-Ny)/Ly
kz = 2.0d0*pi*float((k-1)-Nz)/Lz
phik(i,j,k) = rhok(i,j,k)/( kx*kx + ky*ky + kz*kz)
phik_dum(i,j,k) = phik(i,j,k)
enddo
enddo
enddo
call dfftw_plan_dft_c2r_3d_ (plan_backward, Nx, Ny, Nz, 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
do k = 1, Nz
phi(i,j,k) = phi(i,j,k)/(float(Nx)*float(Ny)*float(Nz))
write(10,*) x(i),y(j),z(k),rho(i,j,k),phi(i,j,k)
enddo
end do
end do
end program poisson3d