-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathrun_grid_test.py
More file actions
157 lines (129 loc) · 5.89 KB
/
run_grid_test.py
File metadata and controls
157 lines (129 loc) · 5.89 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
'''
This program runs hyperion on a grid or predefined models
'''
import itertools
from astropy.table import Table
import models as m
import pickle,os
import numpy as np
from hyperion.util.constants import rsun, lsun, au, msun, yr, c, pc, sigma, G
### these are the parameters to probe:
disk = ["Flared"]#["Alpha"]
name = ["flared_disk_Opt1600"]#['alpha_disk']
folder = ["flared_disk/"]#['alpha_disk/']
T_list=[6000.]
m_sun_list = [0.7]
l_sun_list = [10.]
R_star = np.sqrt(l_sun_list[0]*lsun/(4.0*np.pi*sigma*T_list[0]**4))
disk_mass_list = [0.01]
disk_rmax_list = [100.]
disk_rmin_list = [1]
disk_h0_list = [0.1]
mdot = 2.*lsun*(disk_rmin_list[0]*au/(G*m_sun_list[0]*msun))/3.*2.
print R_star/rsun
#mdot=1e-6
print "mdot = ",mdot/msun*yr
env_mass_list = [1.]
env_rmax_list = [10000.]
env_rmin_list = [50.]
env_power_list = [-2.0]
env_type_list = ['power']
env_rc_list = [100.]
mdot_list = [mdot]
cav_theta_list = [25.]
innerdust_list = ['OH5.hdf5']
outerdust_list = ['OH5.hdf5']
# this creates the list of lists
list_of_lists = [name, folder, l_sun_list, m_sun_list, env_mass_list, env_type_list,disk_rmax_list, env_rmin_list,env_rc_list,env_power_list, disk_h0_list, disk_mass_list , cav_theta_list, disk_rmin_list ,mdot_list, env_rmax_list, T_list,innerdust_list,outerdust_list,disk]
# flatten the list of lists
list_of_lists_flattened = list(itertools.product(*list_of_lists))
# if file exists, load existing grid file
filename = folder[0]+name[0]+'.grid.dat'
if os.path.exists(filename):
grid = pickle.load(open(filename,'r'))
# load up last line of grid
Ngrid = int(grid['name'][-1].split('_')[-1].split('.')[0])
print Ngrid
num = 0
L=len(grid)
for mylist in list_of_lists_flattened:
add = 'yes'
for n in range(L):
curlist = [grid[n][col] for col in ['folder','L_sun','M_sun','env_mass','env_type','disk_rmax','env_rmin','rc','env_power','disk_h0','disk_mass',\
'cav_theta','disk_rmin','mdot','env_rmax','T','innerdustfile','outerdustfile','disk']]
print "Curlist is:",set(mylist[1:])
print "Mylist is:",set(curlist)
print set(mylist[1:]).intersection(set(curlist))
if len(set(mylist[1:]).intersection(set(curlist))) == len(curlist):
add = 'no'
# only add the new row when it's not already in there
if add =='yes':
print "Adding row"
print mylist
grid.add_row(mylist)
num+=1
print "Ngrid+num",Ngrid+num
# change all the names of the new lines
grid['name'][Ngrid+num] = grid['name'][Ngrid+num]+"_"+str(Ngrid+num)
print grid['name'][Ngrid+num]
# save new grid before doing simulation
grid.write(folder[0]+name[0]+".grid.txt",format='ascii.fixed_width')
pickle.dump(grid,open(folder[0]+name[0]+".grid.dat",'wb'))
if num>0:
for i in range(Ngrid+num,len(grid)):
amb_dens = 1e-20
if grid['innerdustfile'][i] == 'd03_5.5_3.0_A.hdf5':
grid['disk_mass'][i] /= 100.
if grid['outerdustfile'][i] == 'd03_5.5_3.0_A.hdf5':
grid['env_mass'][i] /= 100.
amb_dens /= 100.
#if grid['grid['env_type'][i]'][i]==1:
# qdisk = True
#else:
# qdisk=False
model = m.YSOModelSim(name=grid['name'][i],folder=grid['folder'][i],L_sun=grid['L_sun'][i],M_sun=grid['M_sun'][i],env_mass=grid['env_mass'][i],\
env_type=grid['env_type'][i],disk_rmax=grid['disk_rmax'][i],env_rmin=grid['env_rmin'][i],rc=grid['rc'][i],env_power=grid['env_power'][i],\
disk_h_0=grid['disk_h0'][i],disk_mass=grid['disk_mass'][i],cav=True,cav_theta=grid['cav_theta'][i],disk_rmin=grid['disk_rmin'][i],\
mdot=grid['mdot'][i],amb_rmin=1,amb_rmax=grid['env_rmax'][i],env_rmax=grid['env_rmax'][i],T=grid['T'][i],innerdustfile=grid['innerdustfile'][i],\
amb_dens=amb_dens,angles=[0.,10.,20.,30.,40.,50.,60.,70.,80.,90.],angles2=[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],env=True,disk=grid['disk'][i],\
outerdustfile=grid['outerdustfile'][i])
model.initModel()
model.runModel()
model.plotSim(extinction=0,show=True,inc=0,dist_pc=100)
# this just applies to the first time we do a grid
else:
# this creates a grid in the form of a table
allparams = ['name', 'folder','L_sun','M_sun','env_mass','env_type','disk_rmax','env_rmin','rc','env_power','disk_h0','disk_mass',\
'cav_theta','disk_rmin','mdot','env_rmax','T','innerdustfile','outerdustfile','disk']
types = ['<S30','<S30','f8','f8','f8','<S30','f8','f8','f8','f8','f8','f8','f8','f8','f8','f8','f8','<S30','<S30','<S30']
grid = Table(names=allparams,dtype=types)
# only run things if the simulation has not yet been done
for mylist in list_of_lists_flattened:
grid.add_row(mylist)
print "There will be ",len(grid),"models that will be run"
for j in range(len(grid)):
grid['name'][j] = grid['name'][j]+"_"+str(j)
print grid
grid.write(folder[0]+name[0]+".grid.txt",format='ascii.fixed_width')
pickle.dump(grid,open(folder[0]+name[0]+".grid.dat",'wb'))
for i in range(len(grid)):
amb_dens = 1e-20
if grid['innerdustfile'][i] == 'd03_5.5_3.0_A.hdf5':
grid['disk_mass'][i] /= 100.
if grid['outerdustfile'][i] == 'd03_5.5_3.0_A.hdf5':
grid['env_mass'][i] /= 100.
amb_dens /= 100.
#if grid['grid['env_type'][i]'][i]==1:
# qdisk = True
#else:
# qdisk=False
model = m.YSOModelSim(name=grid['name'][i],folder=grid['folder'][i],L_sun=grid['L_sun'][i],M_sun=grid['M_sun'][i],env_mass=grid['env_mass'][i],\
env_type=grid['env_type'][i],disk_rmax=grid['disk_rmax'][i],env_rmin=grid['env_rmin'][i],rc=grid['rc'][i],env_power=grid['env_power'][i],\
disk_h_0=grid['disk_h0'][i],disk_mass=grid['disk_mass'][i],cav=True,cav_theta=grid['cav_theta'][i],disk_rmin=grid['disk_rmin'][i],\
mdot=grid['mdot'][i],amb_rmin=1,amb_rmax=grid['env_rmax'][i],env_rmax=grid['env_rmax'][i],T=grid['T'][i],innerdustfile=grid['innerdustfile'][i],\
amb_dens=amb_dens,angles=[0.,10.,20.,30.,40.,50.,60.,70.,80.,90.],angles2=[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],env=True,disk=grid['disk'][i],\
outerdustfile=grid['outerdustfile'][i])
model.initModel()
model.runModel()
#print grid[i]
model.plotSim(extinction=0,show=True,inc=0,dist_pc=100)