Commit: 8f80def58fb25e3c117d63fa251bac5fe5fa4db8 - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 8f80def58fb25e3c117d63fa251bac5fe5fa4db8 (tree)
• Zeit: 2009-01-20 03:05:19
• Autor: iselllo
• Commiter: iselllo

Log Message

Now the code reads most of the parameters from an input data file (apart from the monomers radia).
The code is now able to save also some intermediate results during the evolution of the system.
Still being developed, check everything!

Ändern Zusammenfassung

• modified: Python-codes/wiener.py

Diff

diff -r dd33bff1ad74 -r 8f80def58fb2 Python-codes/wiener.py
--- a/Python-codes/wiener.py	Mon Jan 19 14:55:24 2009 +0000
+++ b/Python-codes/wiener.py	Mon Jan 19 18:05:19 2009 +0000
@@ -3,11 +3,56 @@
 import numpy as n
 import pylab as p
 
+import csv
+
+
 
 #import absorbe as abso  #compiled fortran code for detecting molecule-to-gas collisions
 import absorbe_snap as absnap  #compiled fortran code for detecting molecule-to-gas collisions
 
 
+#! /usr/bin/env python
+import scipy as s
+import csv
+
+
+def read_input(filename):
+
+
+    param_name_list=[]
+    param_value_list=[]
+
+
+
+    reader = csv.reader(open(filename, 'r'), delimiter = ',')
+    line = 0
+    # I changed the two following lines into  comments since I do not have a header I should skip
+    
+    #IndexHeader = reader.next()
+    #line = line+1
+    for line in reader:
+        #print line
+        param_name, param_value = line[0], line[1]
+
+        param_name_list.append(param_name)
+        param_value_list.append(float(param_value))
+
+
+
+    #print "param_name_list is, ", param_name_list
+    #print "param_name_list is, ", param_value_list
+
+    param_arr=s.asarray(param_value_list)
+
+    #print "param_arr is, ", param_arr
+
+    #print "param_name_list[2] is, ", param_name_list[2]
+
+    return param_arr
+
+
+
+
 
 def evolve_by_time_step(amplitude,t_step,N_part,L,pos_at_t):
 
@@ -39,7 +84,9 @@
 
     n_iter=int(round(t_fin/t_step))
 
-    my_mod=int(round(t_fin/t_step))/10
+    #my_mod=int(round(t_fin/t_step))/10
+
+    save_every=n_iter/number_save
 
     #print "N_iter is, ", n_iter
 
@@ -48,7 +95,7 @@
     
     for i in xrange(n_iter):
 
-        if (s.remainder(i,my_mod) == 0 ):
+        if ((s.remainder(i,save_every) == 0) & (i>0) ):
         
             print "i is, ", i
         #print "monomer_pos are, ", monomer_pos[:,0]
@@ -111,6 +158,149 @@
 
     return results
 
+def evolution_and_absorption_intermediate_save(monomer_pos,number_monomers,amplitude,t_step,\
+                             N_part,L, pos_at_t,r1_arr,ini_coll_mat, number_save):
+
+
+    n_iter=int(round(t_fin/t_step))
+
+
+
+    save_every=n_iter/number_save
+
+    snap_num=0
+
+    
+    for i in xrange(n_iter):
+
+  
+        if (i==0):
+
+            results=absnap.absorption_snap_calc(pos_at_t, monomer_pos[:,0], \
+                                                monomer_pos[:,1], \
+                                                monomer_pos[:,2], r1_arr, i, ini_coll_mat)
+        
+            results_temp=s.copy(results)
+
+
+
+
+        elif (i !=0):
+
+            
+
+            temp=1 #which means I am now saving temporary files during the evolution
+
+
+
+            results=absnap.absorption_snap_calc(pos_at_t, monomer_pos[:,0], \
+                                                monomer_pos[:,1], \
+                                                monomer_pos[:,2], r1_arr, i, results_temp)
+        
+            results_temp=s.copy(results)
+
+            if ((s.remainder(i,save_every) == 0) & (i>0) ):
+
+                snap_num=snap_num + 1
+        
+                print "saving snapshot number , ", snap_num
+
+                #Now save the intermediate results
+
+                abso_temp=absorption_individual_monomers(number_mon,results,time, temp, snap_num )
+
+            
+                stat_coll=stat_collisions(results, time, 0)
+
+                prefix="temp_%01d"%snap_num
+                
+                filename="_collision_statistics.dat"
+
+                filename=prefix+filename
+                
+                p.save(filename,stat_coll )
+
+
+
+
+
+                #Now we perform some statistical analysis of the data
+
+                number_molecules_vs_time=count_molecules(results,N_part,stat_coll, t_step)
+
+                filename="_number_molecules_vs_time.dat"
+
+                filename=prefix+filename
+
+                p.save(filename,number_molecules_vs_time)
+
+                p.save("intermediate_mol_position.dat",pos_at_t )
+                p.save("intermediate_time.dat",(t_step*i*s.ones(1)) )
+
+
+
+        pos_at_t=evolve_by_time_step(amplitude,t_step,N_part,L,pos_at_t)
+
+
+
+
+    return results
+
+
+
+
+
+def absorption_individual_monomers(number_mon,new_abso,time, temp, snap_num):
+
+    for m in xrange(number_mon):
+
+        abso_mon=select_absorption(new_abso,(m+1))  #careful about the m+1 appearing in several places:
+        #it comes from
+        # the fortran routine way of labelling the molecules 
+
+        if (temp ==1):
+
+            prefix="temp_%01d"%snap_num
+                            
+            filename="_abso_%01d"%(m+1)
+
+            filename=prefix+filename
+
+        elif (temp != 1):
+
+            filename="abso_%01d"%(m+1)
+            
+        filename=filename+".dat"
+
+        p.save(filename,abso_mon )
+
+
+
+
+
+        stat_coll=stat_collisions(abso_mon, time, 0)
+
+        if (temp == 1):
+
+            prefix="temp_%01d"%snap_num
+
+            filename="_collision_statistics_%01d"%(m+1)
+
+            filename=prefix+filename
+
+
+        elif (temp !=1):
+
+            filename="collision_statistics_%01d"%(m+1)
+
+
+        filename=filename+".dat"
+
+        p.save(filename,stat_coll )
+
+        return 0
+
+
 
 
 
@@ -481,24 +671,59 @@
 #######################################################################
 #######################################################################
 
-L=3.
-
-t_fin=20.
+# L=3.
 
-t_step=0.001
+# t_fin=20.
 
-amplitude=2.
+# t_step=0.001
 
-N_part=30000
+# amplitude=2.
+
+# N_part=30000
 
 
-number_mon=2  #indicate the number of monomers
+# number_mon=2  #indicate the number of monomers
 
-mon_rad=s.ones(number_mon)*0.5
+# mon_rad=s.ones(number_mon)*0.5
 
 
 
-fold=1  #whether I should fold the trajectories inside the box or they have already been saved that way
+input_arr=read_input("input.csv")
+
+#print "input_arr is, ", input_arr
+
+
+L=input_arr[0]
+
+#print "L is, ", L
+
+t_ini=input_arr[1]
+
+t_fin=input_arr[2]
+
+t_step=input_arr[3]
+
+
+
+N_part=int(round(input_arr[4]))
+number_mon=int(round(input_arr[5]))
+number_save=int(round(input_arr[6]))  #how many intermediate results I am supposed to save
+
+read_ini_pos=int(round(input_arr[7])) #whether I should read the initial position of the molecules \
+# from another simulation or not
+
+read_radia=int(round(input_arr[8])) #whether I should read the list of the radia of the monomers from a file
+
+if (read_radia == 1):
+    mon_rad=p.load("radia.csv", delimiter=",")
+
+
+amplitude=input_arr[9]
+
+
+#I do not need that any longer since I am now using the code in a different way
+
+#fold=1  #whether I should fold the trajectories inside the box or they have already been saved that way
 
 
 
@@ -510,9 +735,9 @@
 time=s.linspace(0.,t_fin,(int(t_fin/t_step)+1) )
 p.save("time.dat", time)
 
+#I am no longe using the "read" parameter: the code no longer reads the trajectories to start with
 
-
-read=1  #it means that the code will read some already-generated trajectories
+#read=1  #it means that the code will read some already-generated trajectories
 # rather than calculate them on  the fly
 
 
@@ -590,8 +815,12 @@
 
 
 
-new_abso=evolution_and_absorption(monomer_pos,number_mon,amplitude,t_step,\
-                             N_part,L, initial_pos,mon_rad ,coll_matrix_ini)
+#new_abso=evolution_and_absorption(monomer_pos,number_mon,amplitude,t_step,\
+#                            N_part,L, initial_pos,mon_rad ,coll_matrix_ini)
+
+
+new_abso=evolution_and_absorption_intermediate_save(monomer_pos,number_mon,amplitude,t_step,\
+                             N_part,L, initial_pos,mon_rad,coll_matrix_ini, number_save)
 
 
 print "OK after test2"