Commit: 22cc1e99ceb6fd3dc01ac219fb45f744a359588d - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 22cc1e99ceb6fd3dc01ac219fb45f744a359588d (tree)
• Zeit: 2008-08-12 00:48:43
• Autor: iselllo
• Commiter: iselllo

Log Message

I fixed a bug in the expression of the kernel (the mass of an aggregate
is NOT equal to the volume, but rather to the number of monomers it
consists of since I am working in units where m_1=1).
FIX THIS BUG also in the code for smoluchowski equation used to evolve
an initially monodisperse aerosol.
I also now calculate the kernel both with <n_i><n_j> and with <n_in_j>.

Ändern Zusammenfassung

• modified: Python-codes/plot_kernel_elements_multiple_directories.py

Diff

diff -r 001ded6119d4 -r 22cc1e99ceb6 Python-codes/plot_kernel_elements_multiple_directories.py
--- a/Python-codes/plot_kernel_elements_multiple_directories.py	Wed Aug 06 18:18:03 2008 +0000
+++ b/Python-codes/plot_kernel_elements_multiple_directories.py	Mon Aug 11 15:48:43 2008 +0000
@@ -6,7 +6,7 @@
 
 time_red=p.load("../1/time_red.dat")  #time I used when comparing the snapshots
 
-n_dir=8  #number of different directories where I will be reading the data from
+n_dir=10  #number of different directories where I will be reading the data from
 
 dir_list=s.arange(n_dir)+1
 
@@ -20,11 +20,11 @@
 print "the snapshot time series is, ", time_red
 
 
-N_mon=5000. #initial number of monomers
+N_mon_tot=5000. #initial number of monomers
 
 rho=0.01
 
-box_vol=N_mon/rho
+box_vol=N_mon_tot/rho
 
 
 by=2
@@ -44,14 +44,21 @@
 
 n_j=s.zeros((n_dir,len(my_selection)))
 
+######################
 
+n_in_j_nonlinear=s.zeros((n_dir,len(my_selection)))    # this array and the one below will be used to calculate <n_in_j> which is NOT <n_i><n_j>.
+
+
+n_in_j_nonlinear_aver=s.zeros(len(my_selection))
+
+################
 n_i_aver=s.zeros(len(my_selection))
 n_j_aver=s.zeros(len(my_selection))
 
 
 sel_ker=s.arange(2)
-sel_ker[0]=1   #NB: these are ORDERED collisions: I can evaluate k_ij for i<j ONLY!!!!
-sel_ker[1]=20
+sel_ker[0]=3   #NB: these are ORDERED collisions: I can evaluate k_ij for i<j ONLY!!!!
+sel_ker[1]=4
 
 
 sel_ker=sel_ker*1.  #to get a floating point array
@@ -134,6 +141,7 @@
 
             n_i[m,i]=len(s.where(size_dist==sel_ker[0])[0])
             n_j[m,i]=len(s.where(size_dist==sel_ker[1])[0])
+            n_in_j_nonlinear[m,i]=n_i[m,i]*n_j[m,i]
 
             #NB: the definition of n_i and n_j is well-posed even if NO collision involving them took place.
 
@@ -146,6 +154,8 @@
             cluster_name=cluster_name+extension
             cluster_name=dir_name+cluster_name
 
+            #print "cluster_name is, ", cluster_name
+            
             recorded_collisions=p.load(cluster_name)
 
 
@@ -219,6 +229,9 @@
 n_i_aver=n_i.mean(axis=0)  #NO ambiguity about these averages; no problem if either n_i or n_j is zero.
 n_j_aver=n_j.mean(axis=0)
 
+n_in_j_nonlinear_aver=n_in_j_nonlinear.mean(axis=0)
+
+n_in_j_nonlinear_aver=n_in_j_nonlinear_aver/box_vol/box_vol #!!!!! careful with the normalization!!
 
 
 number_coll_aver=number_coll.mean(axis=0)
@@ -413,14 +426,71 @@
 
 print "coll_dens is, ", coll_dens
 
+
+
+######## VERY Careful!!!! n_i_nj is now defined as <n_i><n_j> NOT as <n_in_j>.
+
 n_in_j=n_i_aver/box_vol*n_j_aver/box_vol
 
 print "n_in_j is, ", n_in_j
 
 fitted_kernel=residuals_eval(coll_dens, n_in_j, k_ij_inf, k_ij_sup, n_steps)
 
-print "fitted_kernel is, ", fitted_kernel[0]
+if (sel_ker[0]!=sel_ker[1]):
 
+    result=list(fitted_kernel[0])
+    result=s.asarray(result)
+    result=result/2.
+elif (sel_ker[0]==sel_ker[1]):
+    
+    result=list(fitted_kernel[0])
+    result=s.asarray(result)
+#    result=result/2.
+
+    
+
+print "fitted_kernel is, ", result
+
+#print "fitted_kernel is, ", fitted_kernel[0]
+
+
+fitted_kernel_nonlinear=residuals_eval(coll_dens, n_in_j_nonlinear_aver, k_ij_inf, k_ij_sup, n_steps)
+
+if (sel_ker[0]!=sel_ker[1]):
+    result=list(fitted_kernel_nonlinear[0])
+    result=s.asarray(result)
+
+    result=result/2.
+
+elif (sel_ker[0]==sel_ker[1]):
+    
+    result=list(fitted_kernel_nonlinear[0])
+    result=s.asarray(result)
+
+
+print "fitted_kernel is [nonlinear average], ", result
+
+p.save("nonlin_aver_n_i_nj.dat",n_in_j_nonlinear_aver )
+p.save("lin_aver_n_i_nj.dat",n_in_j )
+
+
+
+
+fig = p.figure()
+axes = fig.gca()
+
+
+axes.plot(time_coll,n_in_j, "bo", label="<n_i><n_j>")
+axes.plot(time_coll, n_in_j_nonlinear_aver, "k^", label="<n_in_j>")
+p.xlabel('Time')
+p.ylabel('Number of collisions')
+#p.title("Evolution Mean-free path")
+p.grid(True)
+cluster_name="non_linear_vs_linear_n_in_j.pdf"
+axes.legend()
+p.savefig(cluster_name)
+ 
+p.clf()
 
 
 
@@ -463,6 +533,18 @@
  
 p.clf()
 
+
+
+
+
+p.save("time_coll.dat",time_coll)
+
+p.save("k_ij.dat", fitted_kernel[0])
+
+p.save("coll_ij_density.dat", coll_dens)
+
+p.save("fitted_coll_density.dat",fitted_kernel[0]*n_in_j)
+
 #Now I do as above but I specify a time-span
 
 time_inf=0.
@@ -473,9 +555,21 @@
 
 fitted_kernel=residuals_eval(coll_dens[time_span], n_in_j[time_span], k_ij_inf, k_ij_sup, n_steps)
 
-print "fitted_kernel is, ", fitted_kernel[0]
+    
+if (sel_ker[0]!=sel_ker[1]):
+
+    result=list(fitted_kernel[0])
+    result=s.asarray(result)
+
+    result=result/2.
+
+elif (sel_ker[0]==sel_ker[1]):
+    
+    result=list(fitted_kernel[0])
+    result=s.asarray(result)
 
 
+print "fitted_kernel [on time span] is, ", result
 
 
 fig = p.figure()
@@ -517,4 +611,260 @@
 
 p.save("residuals_span.dat", fitted_kernel[1])
 
+############################################################################################################################
+############################################################################################################################
+############################################################################################################################
+############################################################################################################################
+
+#now I compare the results between the calculated kernel elements, the smoluchowski kernel and "our" kernel
+
+#First some parameters of the system
+
+
+
+beta=1. #cluster/monomer 1/tau
+k_B=1. #in these units
+T_0=0.5 #temperature of the system
+m_mon=1. #monomer mass in these units
+sigma=1. #monomer diameter
+mu=(m_mon*beta)/(3.*s.pi*sigma) # fluid viscosity
+
+
+
+r_mon=0.5
+
+v_mono=4./3.*s.pi*r_mon**3.
+
+print "the monodisperse volume is, ", v_mono
+
+## x=s.logspace(s.log10(v_mono), s.log10(v_mono*4500.),400)   # volume grid
+
+## n_mon=x/v_mono #volume of solids expressed in terms of number of monomers
+
+
+
+
+
+
+
+n_mon=sel_ker*1.  #to have a floating point array.
+
+m=n_mon  # since I set the mass of a single monomer equal to 1, then the mass of a cluster is the same as the number of monomers it consists of.
+
+print "n_mon is, ", n_mon
+
+
+x=n_mon*v_mono
+
+
+threshold=15. # I define the threshold between small and large clusters, to be used where it matters.
+
+
+small=s.where(n_mon<=threshold)
+large=s.where(n_mon>threshold)
+
+
+
+
+a_small=0.385
+df_small=2.25
+
+a_large=0.223
+df_large=1.57
+
+
+a_tot=0.246
+
+df_tot=1.63
+
+
+
+        
+def Brow_ker_cont_optim(Vlist):
+	kern_mat=2.*k_B*T_0/(3.*mu)*(Vlist[:,s.newaxis]**(1./df_tot)+\
+	Vlist[s.newaxis,:]**(1./df_tot))**2./ \
+	(Vlist[:,s.newaxis]**(1./df_tot)*Vlist[s.newaxis,:]**(1./df_tot))
+	return kern_mat
+
+
+def Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs(Vlist):
+	#monomer volume
+	#r_mon=0.5 #monomer radius
+	#v_mon=4./3.*s.pi*r_mon**3.
+	#v_mono already defined as a global variable 
+	#n_mon=Vlist/v_mono #number of monomers in each aggregate
+
+	#print "n_mon is, ", n_mon
+	
+	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
+	#which just depends on the cluster size.
+
+	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
+	#radia of gyration
+
+	#threshold=15.
+	
+# 	small=s.where(n_mon<=threshold)
+# 	large=s.where(n_mon>threshold)
+
+# 	a_small=0.36
+# 	df_small=2.19
+
+# 	a_large=0.241
+# 	df_large=1.62
+
+
+	R_list[small]=a_small*n_mon[small]**(1./df_small)
+
+	R_list[large]=a_large*n_mon[large]**(1./df_large)
+
+#	R_list[0]=0.5   #special case for the monomer radius
+
+#	m=rho_p*Vlist # this holds in general (i.e. for Vlist !=3.)
+	## as long as Vlist is the volume of solid
+	##and not the space occupied by the agglomerate structure
+
+
+
+	c=(8.*k_B*T_0/(s.pi*m))**0.5
+	#print 'c is', c
+	l=8.*Diff/(s.pi*c)
+	#print 'l is', l
+	diam_seq=2.*R_list
+	
+	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
+	
+	beta_fuchs=(\
+	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
+	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
+	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
+	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
+	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
+	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
+	)**(-1.)
+	
+	
+	## now I have all the bits for the kernel matrix
+#	kern_mat=Brow_ker_cont_optim(Vlist)*beta
+
+
+	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
+		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])*beta_fuchs
+
+
+	return kern_mat
+
+
+
+def Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs_explicit_smallest(Vlist):
+	#monomer volume
+	#r_mon=0.5 #monomer radius
+	#v_mon=4./3.*s.pi*r_mon**3.
+	#v_mono already defined as a global variable 
+	#n_mon=Vlist/v_mono #number of monomers in each aggregate
+
+	#print "n_mon is, ", n_mon
+	
+	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
+	#which just depends on the cluster size.
+
+	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
+	#radia of gyration
+
+	#threshold=15.
+	
+# 	small=s.where(n_mon<=threshold)
+# 	large=s.where(n_mon>threshold)
+
+# 	a_small=0.36
+# 	df_small=2.19
+
+# 	a_large=0.241
+# 	df_large=1.62
+
+
+	R_list[small]=a_small*n_mon[small]**(1./df_small)
+
+	R_list[large]=a_large*n_mon[large]**(1./df_large)
+
+        #Now I introduce a modification for R_list < 5
+        
+        sel=s.where(n_mon==1.)
+
+        R_list[sel]=3.87e-1
+
+
+        sel=s.where(n_mon==2.)
+
+        R_list[sel]=6.39e-1
+
+
+        sel=s.where(n_mon==3.)
+
+        R_list[sel]=7.18e-1
+
+        sel=s.where(n_mon==4.)
+
+        R_list[sel]=7.48e-1
+
+
+#	R_list[0]=0.5   #special case for the monomer radius
+
+#	m=rho_p*Vlist # this holds in general (i.e. for Vlist !=3.)
+	## as long as Vlist is the volume of solid
+	##and not the space occupied by the agglomerate structure
+
+       # m=Vlist #since rho = 1.
+
+	c=(8.*k_B*T_0/(s.pi*m))**0.5
+	#print 'c is', c
+	l=8.*Diff/(s.pi*c)
+	#print 'l is', l
+	diam_seq=2.*R_list
+	
+	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
+	
+	beta_fuchs=(\
+	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
+	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
+	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
+	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
+	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
+	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
+	)**(-1.)
+	
+
+        print "Fuchs beta is, ", beta_fuchs
+	
+	## now I have all the bits for the kernel matrix
+#	kern_mat=Brow_ker_cont_optim(Vlist)*beta
+
+
+
+	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
+		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])*beta_fuchs
+
+
+	return kern_mat
+
+
+
+
+smolu_kernel=Brow_ker_cont_optim(x)
+
+print "smoluchowsky kernel is, ", smolu_kernel[0,1]
+
+
+fuchs_kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs(x)
+
+print "fuchs kernel is, ", fuchs_kernel[0,1]
+	
+
+fuchs_modified_below_5=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs_explicit_smallest(x)
+
+print "fuchs_modified_below_5, ", fuchs_modified_below_5[0,1]
+
+
 print "So far so good"
+
+