Commit 15f3220b authored by LAMBERT Jean-charles's avatar LAMBERT Jean-charles
Browse files

remove extra copy/paste

parent bd27ecd6
......@@ -293,7 +293,7 @@ class CCod:
#!!!!sys.exit()
if not self.__fastcod:
#c=CFalcon() # new falcon object
ok,rho,hsml=c.getDensity(pos,mass) # compute density
ok,rho,hsml=CFalcon().getDensity(pos,mass) # compute density
# compute cod
......@@ -369,298 +369,6 @@ class CCod:
return True,comp,i-1
return False,None,None
#
# isMultipleHalo
#
def isMultipleHalo(self):
"""
return if simulation has multiple halo
"""
return self.__is_multiple_halo
#
# ----
#
# compute COD according "select"ion component
def compute(self,select,ncores=None,cod_file=None, fastcod=True,threshold=10000):
"""
compute Center of Density on list of list of components
select: string
list of components, example :
simple list : select=gas,disk
two list : select=gas,disk:stars (note ':' to separate lists)
ncores : integer
#cores used for multiprocessing, if None then all availables cores will be used
cod_file : string
cod's file name, if None it will be automatically build to be stored in sim's directory
fastcod : boolean (default true)
True : compute fas_cod using octree
rue : compute fas_cod using octree
threshold : number of max dens particles estimated by ctree depht keept
if > 0 use this value as number of particles
if < 0 use as percentage of particles
default (10000 particles)
"""
if self.__r is None:
print ("Simulation",self.simname," does not exist",file=sys.stderr)
return False
self.__ctree_threshold=threshold
self.__fastcod=fastcod
if cod_file is None:
self.__cod_file_base=self.__r["dir"]+"/ANALYSIS/"+self.__COD_DIR_NAME
# create directory
if not os.path.isdir(self.__cod_file_base):
try:
os.makedirs( self.__cod_file_base)
except:
print ("Cannot create directory <%s>, aborting"%(self.__cod_file_base),file=sys.stderr)
sys.exit()
if self.__vdebug:
print ("COD FILE BASE =",self.__cod_file_base,file=sys.stderr)
# Figure out best selected components
self.__best_select=""
for colon_s in select.split(":"):
for coma_s in colon_s.split(","):
if self.__vdebug:
print (">> :",coma_s,file=sys.stderr)
if self.__best_select.find(coma_s)==-1: # substring not found
if len(self.__best_select)==0:
sep=""
else:
sep=","
self.__best_select += sep+ coma_s
if self.__vdebug:
print("Best select :",self.__best_select,file=sys.stderr)
# get half of list of snapshot to check if all components exist
half_snap=self.__slist[len(self.__slist)/2]
if self.__vdebug:
print("Half snapshot =",half_snap,file=sys.stderr)
#sys.exit()
test_snap=CSnapshot(half_snap, self.__best_select)
ok=test_snap.nextFrame()
# rebuild select string with existing components
self.__new_select=""
for colon_s in select.split(":"):
ok,data=test_snap.getData(colon_s,"pos")
if not ok:
ok,comp,ii=self.__checkSelectCompN(colon_s)
print(ok,comp,ii,file=sys.stderr)
if ok:
if len(self.__new_select)==0:
sep=""
else:
sep=":"
self.__new_select += sep+ colon_s
if self.__vdebug:
print("new select :",self.__new_select,file=sys.stderr)
test_snap.close()
self.__is_multiple_halo = self.__getHaloParticlesNumber() #
self.__startProcess(select,ncores,cod_file)
self.computeMergingTime()
#
# ----
#
# launch parallel process
def __startProcess(self,select,ncores=None,code_file=None):
"""
Start COD computation in parallel on nores
"""
# compute cores
if ncores is None:
ncores=multiprocessing.cpu_count()
if self.__vdebug:
print ("#cores used :",ncores,file=sys.stderr)
# Feed QUEUE with list of snapshots
q = multiprocessing.Queue()
for snap in self.__slist: # loop on list os snsphots
#print "SNAP put :",snap
q.put(snap) # out a snapshot in the list
# start process jobs
processes=[] # list of processes
n=0
# manage signint signal
#original_sigint_handler = signal.signal(signal.SIGINT, signal.SIG_IGN) # IGNOROE before process creation
# lock to control access to file
lock=Lock()
# loop on all #cores and start a process
for p in range(ncores):
p = Process(target=self.__codProcess, args=(q,n,lock,)) # create
p.start() # start
processes.append(p) # append list of process, used for joining
n += 1
# wait all processes to complete
try:
for p in processes:
print ("waiting.. ",p,file=sys.stderr)
p.join()
except KeyboardInterrupt: # allow to interupt all workers with CTRL+C
for p in processes:
print ("Terminating.. ",p,file=sys.stderr)
p.terminate()
p.join()
while not q.empty():
q.get(block=False)
#
# ----
#
def __codProcess(self,queue_list,n,lock):
"""
Compute COD on core
"""
my_snap=self.__findLastSnapshotWithNoCod(queue_list,n,lock)
stop=False
cpt=0
first=True
while (not stop and my_snap is not None): # and not queue_list.empty() ) :
try:
print ("Core [",n,"] waiting...",file=sys.stderr)
if first: # first time my_snap populated with self.__findLastSnapshotWithNoCod()
first=False
else:
my_snap=queue_list.get(True,0.05) # we must use this line
# True means bloc during 0.05 sec,
# if nothing, then get raise Queue.Empty exception
#my_snap=queue_list.get() # do not use this, could block if nothing to get
#time.sleep(0.01)
uns_snap=CSnapshot(my_snap,self.__best_select)
ok=uns_snap.nextFrame("mxvI")
okt,time=uns_snap.getData("time")
# checkCodTime
for select in self.__new_select.split(":"):
mycod_file=self.__cod_file_base+"/"+self.simname+"."+select+".cod"
lock.acquire()
#exist_time=self.__isTimeInCod(mycod_file,time)
exist_time,tcxv=self.getCodFromFile(mycod_file,time,n)
lock.release()
#print ("MYCOD =",mycod_file,time,exist_time)
if not exist_time:
if self.__vdebug:
print (">> computing density on: ", select,file=sys.stderr)
# check select if it is halo_N string
#
self.__ok_halo_N,self.__halo_comp,self.__halo_N=self.__checkSelectCompN(select)
if self.__ok_halo_N: # then it's an halo_X like
select=self.__halo_comp # rewrite halo_N to halo
if self.__ok_halo_N and not self.__is_multiple_halo:# no multiple Halo, no need to compute
continue
ok,pos=uns_snap.getData(select,"pos")
cxv=np.zeros(6)
if ok and pos.size/3 > 32: # more than 32 particles at that time
ok,mass=uns_snap.getData(select,"mass")
ok,vel =uns_snap.getData(select,"vel")
if self.__ok_halo_N and self.__is_multiple_halo : # WE MUST RE ORDER PARTICLES
print (">> select =",select,file=sys.stderr)
ok,id=uns_snap.getData(select,"id")
print ("ok, id.size=",ok,id.size,file=sys.stderr)
ok,mass,pos,vel=self.__extractHalo(id,mass,pos,vel)
print("after extracting halo mass.size:", mass.size,file=sys.stderr)
#!!!!sys.exit()
if not self.__fastcod:
#c=CFalcon() # new falcon object
ok,rho,hsml=c.getDensity(pos,mass) # compute density
# compute cod
cxv=uns_snap.center(pos,vel,rho*mass)
else : # fastcod
ctree=CTree(pos,vel,mass) # instantiate a ctree object
cxv=ctree.fastCod3(self.__ctree_threshold)
print ("TCXV : ", time,cxv,file=sys.stderr)
else: # no data
pass
lock.acquire()
out=open(mycod_file,"a+")
#cxv=np.insert(cxv,0,time)
#out.write(str(cxv).replace("[","").replace("]","")+"\n")
out.write("%e %e %e %e %e %e %e\n"%(time,cxv[0],cxv[1],cxv[2],cxv[3],cxv[4],cxv[5]))
out.close()
lock.release()
cpt+=1
print ("Core [",n,"] got snap : ",my_snap,cpt,file=sys.stderr)
except Queue.Empty:
stop = True # mo more data
if self.__verbose:
print ("Queue.empty execption trapped...",file=sys.stderr)
#except (KeyboardInterrupt, SystemExit):
# print ("Exiting...")
# #lock.release()
# sys.exit()
# terminate()
# stop = True #
print ("Core [",n,"] DONE !",cpt,file=sys.stderr)
#print "Core [",n,"] got snap : ",queue_list.get()
#
# ----
#
def __extractHalo(self,id,mass,pos,vel):
id_sort=np.argsort(id) # sort according to ID order
offset=0
myidsort=np.zeros(1)
if self.__halo_N==1: # first halo
mylast=self.__halo_part[0]
myidsort = id_sort[0:mylast]
else : # second halo
mylast=self.__halo_part[1]
offset = self.__halo_part[0]+mylast
myidsort = id_sort[self.__halo_part[0]:offset]
#print("myidsort : ",myidsort,mylast,self.__halo_part)
pos = np.reshape(pos,(-1,3)) # reshape pos to nbody,3
pos = np.reshape(pos[myidsort],-1) # select ids
vel = np.reshape(vel,(-1,3)) # reshape vel to nbody,3
vel = np.reshape(vel[myidsort],-1) # select ids
return True,mass[myidsort],pos,vel
#
# ----
#
# check in select string component is indexed like halo_1 halo_2 etc..
def __checkSelectCompN(self,ss,max=10):
"""
check string 'ss' given as parameter is indexed like halo_X (X
"""
xx=ss.find("_")
if (xx!=-1):
comp=ss[0:xx]
#print "comp=",comp
for i in range(max):
if (ss==comp+"_%d"%(i)):
return True,comp,i-1
return False,None,None
#
# ----
#
# try to read model_param.txt file in SIMULATION directory and fill up an array with #particles per galaxy/halo
......@@ -1075,7 +783,7 @@ class CCod:
#!!!!sys.exit()
if not self.__fastcod:
#c=CFalcon() # new falcon object
ok,rho,hsml=c.getDensity(pos,mass) # compute density
ok,rho,hsml=CFalcon().getDensity(pos,mass) # compute density
# compute cod
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment