Homec4science
Diffusion pNbody (master) Browse

pNbody/examples/convolution/nbody-convol.py
c3a510ccee83master

Last Change
Raw File

nbody-convol.py
View Options

#!/usr/bin/env python
import Ptools as pt
from pNbody import *
from pNbody import ic
from pNbody.nbodymodule import *
from numpy import fft
import libcon
n = 100000
eps = 1
Rmax = 10
nb = ic.plummer(n,1,1,1,eps,Rmax,ftype='gadget')
nb.pos = nb.pos*array([1,0,0],float32)
#nb.rename("qq.dat")
#nb.write()
softening = 1e-1
# compute density between 0-N-1
N = 512
xmin = -10.
xmax = 10.
shape = (N,)
val = ones(nb.pos.shape).astype(float32)
mass = nb.mass.astype(float32)
x = nb.x()
x = (x-xmin)/(xmax-xmin)
x = x.astype(float32)
mr = mkmap1d(x,mass,val,shape).astype(float)
# set the kernel
e = softening
e2 = e*e
def Kernel(r):
return -1/sqrt(e2+r*r)
# set the grid
K = 2*N
M = zeros(2*K,float)
G = zeros(2*K,float) # for the "manual convolution"
# set M
for i in xrange(N):
if i>=N:
M[i] = 0
else:
M[i]=mr[i]
for i in xrange(0,N):
g = Kernel(i)
G[i] = g # 0 -> N-1
G[i+K] = g # 2N -> 3N-1
G[K-1 -i] = g # 2N-1 -> N
G[2*K-1 -i] = g # 0 -> N-1
# do the convolution
Phi = zeros(2*K,float)
for i in xrange(0,K): # 0->N-1
for j in xrange(0,K): # 0->N-1
if (i-j)<0:
Gij = G[2*K + i-j]
else:
Gij = G[i-j]
Phi[i] = Phi[i] + Gij * M[j]
Phi = Phi * 4*pi
# do the fourrier transform
G2 = G[:K]
M2 = M[:K]
###################################################
# compute potential using fftw
G2ft = fft.fft(G2)
M2ft = fft.fft(M2)
Phi2ft = G2ft*M2ft
Phi2 = fft.ifft(Phi2ft)
Phi2 = Phi2 * 4*pi
###################################################
# compute potential using fftw
Phi5 = libcon.convolution_1d(G2,M2)
Phi5 = Phi5 * 4*pi
###################################################
# compute potential using cuda
Phi6 = libcon.convolution_1dcu(G2,M2)
Phi6 = Phi6 * 4*pi
print "----",min(Phi)
print "----",min(Phi2)
print "----",min(Phi5)
print "----",min(Phi6)
# compute potential using direct summation
dx = (xmax-xmin)/float(N)
x = arange(xmin,xmax,dx)
pos = zeros((len(x),3),float32)
pos[:,0] = x
#Phi3 = nb.Pot(pos.astype(float32),softening)
Phi3 = zeros(len(x),float)
s2 = softening*softening
for i in xrange(len(x)):
r = x[i] - nb.pos[:,0]
r2 = r*r
Phi3[i] = sum(-nb.mass[:]/sqrt( s2 + r2 ))
pos = zeros((len(x),3),float32)
pos[:,0] = x
Phi4 = nb.Pot(pos.astype(float32),softening)
#Phi4 = nb.Pot(pos.astype(float32),softening)
# convert to int
x = (x-xmin)/(xmax-xmin) * N
#pt.plot(x,Phi3)
#pt.plot(x,Phi4,'--')
#pt.show()
#sys.exit()
#################
# plot
#################
pt.subplot(3,1,1)
pt.plot(M)
pt.plot(M2,'.')
pt.subplot(3,1,2)
pt.plot(G)
pt.plot(G2,'.')
pt.subplot(3,1,3)
pt.plot(Phi)
pt.plot(Phi2,'.')
#pt.plot(x,Phi3,'b--')
#pt.plot(x,Phi4,'c.')
pt.plot(Phi5,'^')
pt.plot(Phi5,'o')
pt.show()
c4science · Help