Commit: 1769ee7ee4c3892a55582399c129c4f9aa815e2f - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 1769ee7ee4c3892a55582399c129c4f9aa815e2f (tree)
• Zeit: 2007-05-22 01:22:46
• Autor: iselllo
• Commiter: iselllo

Log Message

I modified the code Python-codes/postprocess2.py. It now goes reading
the output files to postprocess in an upper directory.
Furthermore, I cleaned up its structures so that it does not do the same
job twice when calculating the mean and saving up the results. I
fundamentally added some if-conditions and slightly changed the way in
which the results are saved into text files. I checked against the
previous versions of the code that it returns the same results as in the
past.

Ändern Zusammenfassung

• modified: Python-codes/postprocess2.py

Diff

diff -r 2446f6bb7d26 -r 1769ee7ee4c3 Python-codes/postprocess2.py
--- a/Python-codes/postprocess2.py	Wed May 16 13:52:40 2007 +0000
+++ b/Python-codes/postprocess2.py	Mon May 21 16:22:46 2007 +0000
@@ -13,12 +13,13 @@
 
 def readq(n1,n2,n3,n_snap_ini,n_snap_fin,name,z_pos,theta_pos,readmean):
 	#q3arr=zeros((n1,n2,n3))
-	q3arr=zeros(n1*n2*n3)
-	qmean=zeros(n1*n2*n3)
-	v_rms=zeros(n1*n2*n3)
-	if (readmean==1):
-			name_bis='means_v_z_3D%03d'%n_snap_fin
-			qmean=pylab.load(name_bis)
+	if (readmean !=1): # I calculate the mean velocity
+		q3arr=zeros(n1*n2*n3)
+	elif (readmean ==1): # I read the mean the velocity & calculate the rms
+		qmean=zeros(n1*n2*n3)
+		v_rms=zeros(n1*n2*n3)
+		name_bis='means_v_z_3D%03d'%n_snap_fin
+		qmean=pylab.load(name_bis)
 	#fluctuations=zeros((n2,(n_snap_fin-n_snap_ini+1)))
 	#flu_av=zeros(n2)
 	jcount=0
@@ -26,54 +27,60 @@
 		print 'l is', l
 		fname = name%l # %03d pads the integers with zeros
 		temp_arr=pylab.load(fname) # loading
-		q3arr[:]=temp_arr[:]+q3arr[:]
-		if (readmean==1):
+		if (readmean != 1): # I define the array I need for the mean velocity
+			q3arr[:]=temp_arr[:]+q3arr[:]
+		elif (readmean==1): # I define the array for the rms & I simply read the\ 
+			#mean velocity I had already calculated
 			v_rms[:]=(temp_arr[:]-qmean[:])*(temp_arr[:]-qmean[:])+v_rms[:]
-			
 		jcount=jcount+1
 		#print 'the dimension of q3arr is', shape(q3arr)
 		#now I divide by the number of snapshots in order to have a mean field
-	q3arr=q3arr/float(n_snap_fin-n_snap_ini+1.) #I first take a temporal mean
-	name_bis='means_v_z_3D%03d'%n_snap_fin
-	pylab.save(name_bis,q3arr)
-	v_rms=v_rms/float(n_snap_fin-n_snap_ini+1.) # temporal mean of v_rms
-	q3arr=reshape(q3arr,[n1,n2,n3])
-	v_rms=reshape(v_rms,[n1,n2,n3])
-	my_mean=zeros(n2)
-	my_v_rms=zeros(n2)
-	
-	
-#The following lines are a much faster and cleaner way to deal with the means along theta 
-# and r
-	my_mean=q3arr.mean(axis=0).mean(axis=1)
-	my_v_rms=v_rms.mean(axis=0).mean(axis=1)
-	my_v_rms=sqrt(my_v_rms)
-	#now I work out the average v_z(r) for a specific z (average on a tube cross
-	#section) 
-	mean_cross=zeros(n2)
-	#for i in range(0,n1):
-		#mean_cross[:]= q3arr[i,:,z_pos]+mean_cross[:]
-	#mean_cross=mean_cross/n1
-	mean_cross=q3arr[:,:,z_pos].mean(axis=0)
-	#finally I work out the temporal mean for a specific theta and x
-	mean_theta=q3arr[theta_pos,:,z_pos]
-	
-	#now I glue together the results
-	my_means=zeros((n2,4))
-	my_means[:,0]=my_mean
-	my_means[:,1]=mean_cross
-	my_means[:,2]=mean_theta
-	my_means[:,3]=my_v_rms
-	name_bis='means_v_z%03d'%n_snap_fin
-	pylab.save(name_bis,my_means)
-	
+	if (readmean != 1):
+		#I simply take the temporal mean and save the array along
+		# the whole tube. It will be used afterwards to
+		#evaluate the rms	
+		q3arr=q3arr/float(n_snap_fin-n_snap_ini+1.) #I first take a temporal mean
+		name_bis='means_v_z_3D%03d'%n_snap_fin
+		pylab.save(name_bis,q3arr)
+		# Now I conveniently change the shape of the array 
+		q3arr=reshape(q3arr,[n1,n2,n3])
+		#my_mean=zeros(n2)
+		#Now I take the mean along theta and z to get v_z(r)
+		my_mean=q3arr.mean(axis=0).mean(axis=1)
+		#now I work out the average v_z(r) for a specific 
+		#z (average on a tube cross section)
+		#mean_cross=zeros(n2)
+		#similar procedure as before, but along only an axis only
+		mean_cross=q3arr[:,:,z_pos].mean(axis=0)
+		#finally I work out the temporal mean for a specific theta and x
+		mean_theta=q3arr[theta_pos,:,z_pos]
+		#now I glue together the results
+		my_means=zeros((n2,3))
+		my_means[:,0]=my_mean
+		my_means[:,1]=mean_cross
+		my_means[:,2]=mean_theta
+		#my_means[:,3]=my_v_rms
+		name_bis='means_v_z%03d'%n_snap_fin
+		pylab.save(name_bis,my_means)
+	elif (readmean ==1):
+		v_rms=v_rms/float(n_snap_fin-n_snap_ini+1.) # temporal mean of v_rms
+		v_rms=reshape(v_rms,[n1,n2,n3])
+		#my_v_rms=zeros(n2)
+		#The following lines are a much faster and cleaner way to deal
+		# with the means along theta and r
+		my_v_rms=v_rms.mean(axis=0).mean(axis=1)
+		my_v_rms=sqrt(my_v_rms)
+		name_bis='v_z_rms%03d'%n_snap_fin
+		#my_means=zeros(n2)
+		pylab.save(name_bis,my_v_rms)
+		my_means=my_v_rms
 	return my_means
 
 
 
-n1=65 # theta grid
-n2=65 # r grid
-n3=65 # z grid
+n1=129 # theta grid
+n2=88 # r grid
+n3=129 # z grid
 
 z_pos=25 # I fix a certain position along z
 # I am not really using this now, but it could come handy later on.
@@ -83,7 +90,8 @@
 n_snap_fin = int(raw_input("Enter the value of the final snapshot to read    "))
 #n_snap_ini=15
 #n_snap_fin=20 # number of snapshots
-name='vz%05d'
+name='../vz%05d'
+#print 'name is ', name
 readmean = int(raw_input("Do you want to read (1) qmean for the v_rms?    "))
 
 
@@ -97,7 +105,7 @@
 elif (icalc == 0):
 	print 'here I am'
 	#l=n_snap_fin
-	fname='means_v_z%03d'%n_snap_fin
+	fname='../means_v_z%03d'%n_snap_fin
 	print 'fname is', fname
 	my_arr=pylab.load(fname)
 
@@ -109,7 +117,7 @@
 
 
 
-data=pylab.load("rg2.out")
+data=pylab.load("../rg2.out")
 #data=pylab.load("radstr.out")
 
 print 'data[:,0] is', data[:,0]
@@ -126,19 +134,21 @@
 	vel_prof=(1-my_arr)**(1./6.)
 	return vel_prof
 
-
-name_bis='Average_v_z(r)_on_snapshot%03d'%n_snap_fin
+if (readmean != 1): # in this case I plot some statistics about the average velocity I
+	# have just worked out
 
-vel_emp=one_sixth(data[:,0])
-pylab.plot(data[:,0],my_arr[:,0],data[:,0],my_arr[:,1],data[:,0],my_arr[:,2]\
-,data[:,0],vel_emp*max(my_arr[:,0]))
-pylab.legend(('average on the whole tube','average on cross section','specific position'\
-,'empirical correlation'))
-pylab.xlabel('Radial Coordinate')
-pylab.ylabel('Physical Coordinate')
-pylab.title('Radial Grid')
-pylab.grid(True)
-pylab.savefig(name_bis)
+	name_bis='Average_v_z(r)_on_snapshot%03d'%n_snap_fin
+
+	vel_emp=one_sixth(data[:,0])
+	pylab.plot(data[:,0],my_arr[:,0],data[:,0],my_arr[:,1],data[:,0],my_arr[:,2]\
+	,data[:,0],vel_emp*max(my_arr[:,0]))
+	pylab.legend(('average on the whole tube','average on cross section',\
+	'specific position','empirical correlation'))
+	pylab.xlabel('Radial Coordinate')
+	pylab.ylabel('Physical Coordinate')
+	pylab.title('Radial Grid')
+	pylab.grid(True)
+	pylab.savefig(name_bis)
 
 # now I compare everything with the result I got with Comsol
 
@@ -147,36 +157,38 @@
 #comsol[:,0]=comsol[:,0]/max(comsol[:,0])
 #comsol[:,1]=comsol[:,1]/max(comsol[:,1])*max(my_arr[:,0])
 
-name_bis='Average_v_z(r)_on_snapshot%03d_bis'%n_snap_fin
+	name_bis='Average_v_z(r)_on_snapshot%03d_bis'%n_snap_fin
 
-pylab.plot(data[:,0],my_arr[:,0],\
-data[:,0],vel_emp*max(my_arr[:,0]))
-pylab.legend(('average on the whole tube',\
-'empirical correlation'))
-pylab.xlabel('Radial Coordinate')
-pylab.ylabel('Physical Coordinate')
-pylab.title('Radial Grid')
-pylab.grid(True)
-pylab.savefig(name_bis)
-
+	pylab.plot(data[:,0],my_arr[:,0],\
+	data[:,0],vel_emp*max(my_arr[:,0]))
+	pylab.legend(('average on the whole tube',\
+	'empirical correlation'))
+	pylab.xlabel('Radial Coordinate')
+	pylab.ylabel('Physical Coordinate')
+	pylab.title('Radial Grid')
+	pylab.grid(True)
+	pylab.savefig(name_bis)
 
 
-name_bis='v_z(r)_rms%03d'%n_snap_fin
+if (readmean ==1): # in this case I am ready with the rms plots
+	
 
-pylab.plot(data[:,0],my_arr[:,3])
-#pylab.legend(('average on the whole tube'))
-pylab.xlabel('Radial Coordinate')
-pylab.ylabel('Physical Coordinate')
-pylab.title('Radial Grid')
-pylab.grid(True)
-pylab.savefig(name_bis)
+	name_bis='v_z(r)_rms%03d'%n_snap_fin
 
-name_bis='U_z_rms%03d'%n_snap_fin
+	pylab.plot(data[:,0],my_arr)
+	#pylab.legend(('average on the whole tube'))
+	pylab.xlabel('Radial Coordinate')
+	pylab.ylabel('Physical Coordinate')
+	pylab.title('Radial Grid')
+	pylab.grid(True)
+	pylab.savefig(name_bis)
 
-results=zeros((n2,2))
-results[:,0]=data[:,0]
-results[:,1]=my_arr[:,3]
-pylab.save(name_bis,results)
+	name_bis='U_z_rms%03d'%n_snap_fin
+
+	results=zeros((n2,2))
+	results[:,0]=data[:,0]
+	results[:,1]=my_arr
+	pylab.save(name_bis,results)
 
 #name_bis='Radial_velocity_fluctuations%03d'%n_snap_fin
 
@@ -191,24 +203,23 @@
 
 
 #Now it is time to calculate the physical axial velocity in the tube
-
-U_phys=my_arr[:,0]*U_b*2.
+if (readmean != 1):
+	U_phys=my_arr[:,0]*U_b*2.
 
 
-#print 'U_phys is',U_phys
-
-#r_coord=r_ini
-r_phys=r_ini*D/2.
+	#print 'U_phys is',U_phys
 
-name_bis='physical_axial_velocity%03d'%n_snap_fin
-pylab.plot(r_phys,U_phys)
-#pylab.legend(('average on the whole tube','comsol',\
-#'empirical correlation'))
-pylab.xlabel('Radial Coordinate')
-pylab.ylabel('Physical Coordinate')
-pylab.title('Axial Velocity vs r')
-pylab.grid(True)
-pylab.savefig(name_bis)
+	#r_coord=r_ini
+	r_phys=r_ini*D/2.
+	name_bis='physical_axial_velocity%03d'%n_snap_fin
+	pylab.plot(r_phys,U_phys)
+	#pylab.legend(('average on the whole tube','comsol',\
+	#'empirical correlation'))
+	pylab.xlabel('Radial Coordinate')
+	pylab.ylabel('Physical Coordinate')
+	pylab.title('Axial Velocity vs r')
+	pylab.grid(True)
+	pylab.savefig(name_bis)
 
 
 # Now I work out different units
@@ -234,50 +245,54 @@
 print 'delta_nu is', delta_nu
 
 #Now the same plot as before in wall units (for r and v_z)
-r_plus=r_phys/delta_nu
-U_plus=U_phys/u_tau
-
-name_bis='axial_velocity_wall_units%03d'%n_snap_fin
-pylab.plot(r_plus,U_plus)
-#pylab.legend(('average on the whole tube','comsol',\
-#'empirical correlation'))
-pylab.xlabel('r+')
-pylab.ylabel('u+')
-pylab.title('Axial Velocity y+ vs r+')
-pylab.grid(True)
-pylab.savefig(name_bis)
-
-r_red=r_plus[2:(len(r_ini)-1)]
+if (readmean != 1):
+	
+	r_plus=r_phys/delta_nu
+	U_plus=U_phys/u_tau
 
-name_bis='axial_velocity_wall_units_distance_from_wall%03d'%n_snap_fin
-pylab.plot(max(r_plus)-r_plus,U_plus/max(U_plus),max(r_plus)-r_plus,max(r_plus)-r_plus,\
-max(r_plus)-r_red,\
-(5.2+(1./0.41)*log(max(r_plus)-r_red))/max((5.2+(1./0.41)*log(max(r_plus)-r_red))))
-pylab.axis([0.,max(r_plus),0.,1.])
-pylab.legend(('numerics','linear behavior','log-law'))
-pylab.xlabel('r+')
-pylab.ylabel('u+')
-pylab.title('Axial Velocity y+ vs r+')
-pylab.grid(True)
-pylab.savefig(name_bis)
+	name_bis='axial_velocity_wall_units%03d'%n_snap_fin
+	pylab.plot(r_plus,U_plus)
+	#pylab.legend(('average on the whole tube','comsol',\
+	#'empirical correlation'))
+	pylab.xlabel('r+')
+	pylab.ylabel('u+')
+	pylab.title('Axial Velocity y+ vs r+')
+	pylab.grid(True)
+	pylab.savefig(name_bis)
 
-name_bis='axial_velocity_wall_units_distance_from_wall_unscaled%03d'%n_snap_fin
-pylab.plot(max(r_plus)-r_plus,U_plus,max(r_plus)-r_plus,max(r_plus)-r_plus,\
-max(r_plus)-r_red,\
-(5.2+(1./0.41)*log(max(r_plus)-r_red)))
-pylab.axis([0.,max(r_plus),0.,max(max(U_plus),\
-max((5.2+(1./0.41)*log(max(r_plus)-r_red))))])
-pylab.legend(('numerics','linear behavior','log-law'))
-pylab.xlabel('r+')
-pylab.ylabel('u+')
-pylab.title('Axial Velocity y+ vs r+')
-pylab.grid(True)
-pylab.savefig(name_bis)
+	r_red=r_plus[2:(len(r_ini)-1)]
 
-name_bis='r_and_Uz%03d'%n_snap_fin
-#results=zeros((len(U_plus),2))
-results[:,0]=data[:,0]
-results[:,1]=my_arr[:,0]
-pylab.save(name_bis,results)
+	name_bis='axial_velocity_wall_units_distance_from_wall%03d'%n_snap_fin
+	pylab.plot(max(r_plus)-r_plus,U_plus/max(U_plus),\
+	max(r_plus)-r_plus,max(r_plus)-r_plus,\
+	max(r_plus)-r_red,\
+	(5.2+(1./0.41)*log(max(r_plus)-r_red))/max((5.2+(1./0.41)*\
+	log(max(r_plus)-r_red))))
+	pylab.axis([0.,max(r_plus),0.,1.])
+	pylab.legend(('numerics','linear behavior','log-law'))
+	pylab.xlabel('r+')
+	pylab.ylabel('u+')
+	pylab.title('Axial Velocity y+ vs r+')
+	pylab.grid(True)
+	pylab.savefig(name_bis)
+
+	name_bis='axial_velocity_wall_units_distance_from_wall_unscaled%03d'%n_snap_fin
+	pylab.plot(max(r_plus)-r_plus,U_plus,max(r_plus)-r_plus,max(r_plus)-r_plus,\
+	max(r_plus)-r_red,\
+	(5.2+(1./0.41)*log(max(r_plus)-r_red)))
+	pylab.axis([0.,max(r_plus),0.,max(max(U_plus),\
+	max((5.2+(1./0.41)*log(max(r_plus)-r_red))))])
+	pylab.legend(('numerics','linear behavior','log-law'))
+	pylab.xlabel('r+')
+	pylab.ylabel('u+')
+	pylab.title('Axial Velocity y+ vs r+')
+	pylab.grid(True)
+	pylab.savefig(name_bis)
+
+	name_bis='r_and_Uz%03d'%n_snap_fin
+	results=zeros((len(U_plus),2))
+	results[:,0]=data[:,0]
+	results[:,1]=my_arr[:,0]
+	pylab.save(name_bis,results)
 
 print 'So far so good'
\ No newline at end of file