Page MenuHomec4science
Files F121882488

test_mpi.py
No OneTemporary
Actions

File Metadata

Created
Mon, Jul 14, 14:28

test_mpi.py
View Options

#!/usr/bin/env python
from pNbody import *
from pNbody import ic
from pNbody import libutil
from numpy import *
import sys
import time
import copy
from PyGear import gadget
nb = Nbody("snap.dat",ftype="gadget")
###################################
# init PyGear
gadget.InitMPI() # init MPI
gadget.InitDefaultParameters() # init default parameters
gadget.Info()
params = {}
params['ErrTolTheta'] = 0.7
params['DesNumNgb'] = 15
params['MaxNumNgbDeviation'] = 3
params['UnitLength_in_cm'] = 3.085e+21
params['UnitMass_in_g'] = 4.435693e+44
params['UnitVelocity_in_cm_per_s'] = 97824708.2699
params['BoxSize'] = 1.0
params['PartAllocFactor'] = 10
gadget.SetParameters(params)
###################################
# load particles
gadget.LoadParticles(array(nb.npart),nb.pos,nb.vel,nb.mass,nb.num,nb.tpe)
###################################
# compute Delaunay
t1 = time.time()
gadget.ConstructDelaunay()
t2 = time.time()
print "Delaunay","in",(t2-t1),"s"
###################################
# compute Voronoi
t1 = time.time()
gadget.ComputeVoronoi()
t2 = time.time()
print "Voronoi","in",(t2-t1),"s"
#sys.exit()
###################################
# retrieve ghost points
gpos = gadget.GetAllGhostPositions()
gnb = Nbody(pos=gpos,ftype='gadget')
gnb.u = zeros(gnb.nbody)
gnb.rho = gadget.GetAllGhostvDensities()
gnb.rsp = gadget.GetAllGhostvVolumes()
gnb.init()
# save
gnb.rename('ghost.dat')
gnb.set_pio("yes")
gnb.write()
###################################
# retrieve points
nb.pos = gadget.GetAllPositions()
nb.vel = gadget.GetAllVelocities()
nb.mass = gadget.GetAllMasses()
nb.num = gadget.GetAllID()
nb.tpe = None
nb.u = gadget.GetAllEnergySpec() # not defined
nb.rho = gadget.GetAllvDensities() # note the v
nb.rsp = gadget.GetAllvVolumes()
nb.rsp = gadget.GetAllDensities() # put sph densities in rsp
nb.init()
# save
nb.rename('snap2.dat')
nb.set_pio("yes")
nb.write()
# save all
#nbtot = nb + gnb
nbtot = copy.deepcopy(nb) # this is important to avoid to mix particles
nbtot.append(gnb,do_not_sort=True) # then they are incompatible with vpos below
nbtot.rename('snap+ghost.dat')
nbtot.write()
###################################
# get some values
TriangleList = gadget.GetAllDelaunayTriangles()
rho,mn,mx,cd = libutil.set_ranges(nbtot.rho,scale='log')
rho = rho/255.
#rho = clip(nbtot.rho,0,1)
###################################
# plot
###################################
#sys.exit()
import Ptools as pt
import plot
########################
# draw a box
plot.draw_box(x=array([0,1,1,0]),y=array([0,0,1,1]))
########################
# draw the triangles
i = 0
for Triangle in TriangleList:
P1 = Triangle['coord'][0]
P2 = Triangle['coord'][1]
P3 = Triangle['coord'][2]
#plot.draw_triangle(P1,P2,P3,c='r')
#cm = 1/3.*(P1+P2+P3)
#pt.text(cm[0],cm[1],Triangle['id'],fontsize=12,horizontalalignment='center',verticalalignment='center')
########################
# draw the points
#pt.scatter( nb.pos[:,0], nb.pos[:,1],lw=0,s=50,color='r')
#pt.scatter(gnb.pos[:,0],gnb.pos[:,1],lw=0,s=50)
########################
# draw the Voronoi
cmap = pt.GetColormap('rainbow4',revesed=False)
# single point with vpoints
for i in xrange(nbtot.nbody):
vpos = gadget.GetvPointsForOnePoint(i)
#pt.scatter(nb.pos[i,0],nb.pos[i,1],lw=0,s=50)
plot.draw_cell(vpos,color=[rho[i],rho[i],rho[i]],alpha=0.8)
pt.axis([-0.1,1.1,-0.1,1.1])
pt.show()

Event Timeline

c4science · Help