Peano
Loading...
Searching...
No Matches
para_plot.py
Go to the documentation of this file.
1import numpy as np
2import matplotlib.pyplot as plt
3from scipy import io
4import os
5import time
6import argparse
7import math
8
9
10pi=np.pi
11G=6.67408*1e-11
12M=0.15
13rho_c_ini=1
14rho_c0_value=1.8788*10**(-26)
15a_i=0.001
16
17def file_len(fname):
18 with open(fname) as f:
19 for i, l in enumerate(f):
20 pass
21 return i + 1
22
23#lis=os.listdir("./")
24#for item in lis:
25# if item.endswith(".pvd") or item.endswith(".vtp"): os.remove(os.path.join("./",item))
26
27parser = argparse.ArgumentParser(description='Tracer_conversion')
28parser.add_argument("-if", "--input-file", dest="file_name", required="True", type=str, help="input csv file name from paraview(and directory)" )
29parser.add_argument("-pq", "--plot-quantity", dest="pq", default="rho", type=str, choices=["rho","v_x","v_y","v_z","E_t","P"], help="choose the quantity you want to plot" )
30parser.add_argument("-rescale", "--rescale", dest="IfRescale", action="store_true", default=False, help="if you need the rescaling using turnaround quantites" )
31parser.add_argument("-ylog", "--yscale-log", dest="Ifylog", action="store_true", default=False, help="if you need the plot with y axis in log scale" )
32parser.add_argument("-XUcode", "--x-code-unit", dest="IfXunitInCode", action="store_true", default=False, help="enable this option, if you want the plot in code unit for distance (the default one is in physical unit)" )
33parser.add_argument("-Xrescale", "--x-rescale", dest="IfXRescale", action="store_true", default=False, help="enable this option, if you want the x axis also get rescaled to lambda in Bertschinger, 1985. conflict with XUcode" )
34parser.add_argument("-RawD", "--raw-data", dest="IfRawData", action="store_true", default=False, help="enable this option, if you want the plot of the raw data (in code unit)" )
35args = parser.parse_args()
36
37IfXunitInCode=args.IfXunitInCode
38#file_name=args.file_name
39#dt=args.dt;
40
41
43
44#workaround=args.file_name[5:7]
45#snapshots=[int(workaround)]
46snapshots=[10,30,50]
47snapshot=0
48
49rho = [[]];
50m_1 = [[]]; m_2 = [[]]; m_3 = [[]];
51E_t = [[]]
52coor= [[]]
53for snap in snapshots:
54 f=open("para_"+str(snap)+".csv")
55 print("Read in datafile: "+"para_"+str(snap)+".csv")
56
57 #skip the description line
58 dat=f.readlines()[1:]
59
60 count=0
61 for line in dat:
62 tem=list(map(float,line.split(',')))
63 if math.isnan(tem[1]):
64 continue
65 coor[snapshot].append(tem[7]); rho[snapshot].append(tem[1]);
66 m_1[snapshot].append(tem[2]); m_2[snapshot].append(tem[3]); m_3[snapshot].append(tem[4]);
67 E_t[snapshot].append(tem[5]);
68
69 print("sample points received:", len(rho[snapshot]),"for snapshot code time t=", snap)
70 snapshot+=1
71 rho.append([]);coor.append([])
72 m_1.append([]);m_2.append([]);m_3.append([])
73 E_t.append([])
74 #rho[1].append(12)
75 #print(rho)
76 f.close()
77
78sample_point=len(rho[0])
79origin_point=int(sample_point/2)+1
80
81q_id=args.pq
82
83for snapshot in range(len(snapshots)):
84#==========================================================
85# Read in data
86#==========================================================
87 #t_r_id=tstep-t_id-1
88 #t_code=data_set[0,t_r_id,0]
89 #data_at_t=data_set[:,t_r_id,:]
90 #sort_data=data_at_t[np.argsort(data_at_t[:,1])]
91
92 coor_data=np.array(coor[snapshot])
93 if q_id=="P":
94 e_data = np.array(E_t[snapshot])
95 m1_data = np.array(m_1[snapshot])
96 m2_data = np.array(m_2[snapshot])
97 m3_data = np.array(m_3[snapshot])
98 rho_data = np.array(rho[snapshot])
99 elif q_id=="rho":
100 q_data=np.array(rho[snapshot])
101 elif q_id=="v_x":
102 q_data=np.array(m_1[snapshot])
103 elif q_id=="v_y":
104 q_data=np.array(m_2[snapshot])
105 elif q_id=="v_z":
106 q_data=np.array(m_3[snapshot])
107 elif q_id=="E_t":
108 q_data=np.array(E_t[snapshot])
109 if q_id in ["v_x","v_y","v_z"]:
110 rho_data=np.array(rho[snapshot])
111 #print(sort_data)
112
113#==========================================================
114# start to calculate the real physical quantites
115#==========================================================
116 if not args.IfRawData:
117 t_code=snapshots[snapshot]
118 t_real=(0.1694*a_i**(-0.5)-t_code/11.8)**(-3)*10**15 #*h^(-1)
119 a_scale=0.0287*(0.1694*a_i**(-0.5)-t_code/11.8)**(-2) # when code time ~ 2*(a_i^(-0.5)-1), a~1
120 print("cosmological time:", t_real)
121 print("a_scale:", a_scale)
122
123 coor_data_physical=100*a_scale*coor_data
124
125 if q_id in ["v_x","v_y","v_z"]:
126 q_id_show="$"+q_id+"/(km/s)$"
127 q_data=np.array(q_data)/np.array(rho_data) #remove density from momentum
128 q_data=q_data*10**4/a_scale #H_0*L= 100h*100h^(-1) km/s/Mpc * Mpc = 1e4
129 v_hubble=100*(2.057*10**17/t_real)*coor_data_physical #hubble velocity here
130 q_data=q_data+v_hubble
131 elif q_id=="rho":
132 q_id_show=r"$\rho/\rho_c $"
133 q_data=q_data/(a_scale**3)
134 elif q_id=="P":
135 q_id_show="P/(Pa)"
136 q_data=(5/3-1)*(e_data-0.5*(m1_data**2+m2_data**2+m3_data**2)/rho_data)
137 q_data=q_data*10**(8+6)*rho_c0_value/(a_scale**5)
138 else:
139 q_id_show="$"+q_id+"/(Code Unit)"
140
141#==============================================================
142# now do the rescaling
143#==============================================================
144 r_ta=a_i*100*3.086*1e19*((3*M)/(4*pi*rho_c_ini))**(1/3)*(4/3/pi)**(8/9)*(t_real/6.505/1e12)**(8/9) # km. eq(31), RAMSES
145 if args.IfRescale:
146 if q_id=="rho":
147 q_data=q_data*rho_c0_value # *h^2 kg/m^3
148 rho_scaling=1/(6*pi*G*t_real**2) #h^2 kg/m^3
149 q_data=q_data/rho_scaling #actually it is equal to code density
150 q_id_show=r"$\rho/\rho_{ta}(t) $"
151 if q_id in ["v_x","v_y","v_z"]:
152 #r_ta=a_scale*100*3.086*1e19*((3*M)/(4*pi*rho_c_ini))**(1/3)*(4/3/pi)**(8/9)*(t_real/6.505/1e12)**(8/9) # km. eq(31), RAMSES
153 v_scaling=r_ta/t_real
154 q_data=q_data/v_scaling
155 q_id_show=r"$u/u_{ta}(t) $"
156 print("scaling velocity (km/s):", v_scaling)
157 if q_id=="P":
158 rho_scaling=1/(6*pi*G*t_real**2)
159 P_scaling=rho_scaling*r_ta**2*1e6/(t_real**2)
160 q_data=q_data/P_scaling
161 q_id_show=r"$P/P_{ta}(t) $"
162
163 if IfXunitInCode:
164 plt.plot(coor_data[origin_point:],q_data[origin_point:],label="t= "+str(t_code)+" a= "+str(a_scale))
165 elif args.IfXRescale:
166 plt.plot(np.log10(coor_data_physical[origin_point:]*3.086*1e19/r_ta),q_data[origin_point:],label="t= "+str(t_code)+" a= "+str(a_scale))
167 else:
168 plt.plot(coor_data_physical[origin_point:],q_data[origin_point:],label="t= "+str(t_code)+" a= "+str(a_scale))
169
170#=================================================================
171# or we just plot original raw data
172#=================================================================
173 else:
174 t_code=snapshots[snapshot]
175 a_scale=0.0287*(0.1694*a_i**(-0.5)-t_code/11.8)**(-2) # when code time ~ 2*(a_i^(-0.5)-1), a~1
176 if q_id in ["v_x","v_y","v_z"]:
177 q_data=np.array(q_data)/np.array(rho_data) #remove density from momentum
178 elif q_id=="P":
179 q_data=(5/3-1)*(e_data-0.5*(m1_data**2+m2_data**2+m3_data**2)/rho_data)
180
181 q_id_show="$"+q_id+"/(Code Unit)$"
182 plt.plot(coor_data[origin_point:],q_data[origin_point:],label="t= "+str(t_code)+" a= "+str(a_scale))
183
184plt.title(q_id)
185plt.legend()
186if IfXunitInCode:
187 plt.xlabel("$x/(Code Unit)$")
188elif args.IfXRescale:
189 plt.xlabel(r"$Log_{10}\lambda$")
190else:
191 plt.xlabel("$x/(h^{-1}Mpc)$")
192plt.ylabel(q_id_show)
193#plt.ylim([1.0475,1.0525])
194if args.Ifylog:
195 plt.yscale("log")
196
197#if IfXunitInCode:
198 #plt.xscale("log")
199 #plt.xlim([0.02,0.52])
200
201plt.show()
202
203
204
205
206
207
208
209
210
para_plot.str
str
Definition para_plot.py:28
para_plot.file_len
file_len(fname)
Definition para_plot.py:17