-
Notifications
You must be signed in to change notification settings - Fork 50
/
freefall.py
150 lines (131 loc) · 5.36 KB
/
freefall.py
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
########################################################################
#
# Free-fall example script
#
#
# Copyright (c) 2013-2016, Grackle Development Team.
#
# Distributed under the terms of the Enzo Public Licence.
#
# The full license is in the file LICENSE, distributed with this
# software.
########################################################################
from matplotlib import pyplot
import os
import yt
from pygrackle import \
chemistry_data, \
FluidContainer, \
evolve_constant_density, \
evolve_freefall
from pygrackle.utilities.physical_constants import \
mass_hydrogen_cgs, \
mass_electron_cgs, \
sec_per_Myr, \
cm_per_mpc
tiny_number = 1e-60
if __name__=="__main__":
current_redshift = 0.
# Set solver parameters
my_chemistry = chemistry_data()
my_chemistry.use_grackle = 1
my_chemistry.with_radiative_cooling = 1
my_chemistry.primordial_chemistry = 3
my_chemistry.UVbackground = 0
my_chemistry.self_shielding_method = 0
my_chemistry.H2_self_shielding = 0
my_chemistry.Gamma = 5. / 3.
my_chemistry.CaseBRecombination = 0
my_chemistry.cie_cooling = 1
my_chemistry.h2_optical_depth_approximation = 1
my_chemistry.interstellar_radiation_field = 0.
if os.environ.get("METAL_COOLING", 0) == "1":
my_chemistry.metal_cooling = int(os.environ["METAL_COOLING"])
my_dir = os.path.dirname(os.path.abspath(__file__))
grackle_data_file = bytearray(os.path.join(
my_dir, "..", "..", "..", "input", "cloudy_metals_2008_3D.h5"), 'utf-8')
my_chemistry.grackle_data_file = grackle_data_file
my_chemistry.h2_on_dust = 1
my_chemistry.use_dust_density_field = 1
metallicity = 1e-3
else:
my_chemistry.metal_cooling = 0
# Set units
my_chemistry.comoving_coordinates = 0 # proper units
my_chemistry.a_units = 1.0
my_chemistry.a_value = 1. / (1. + current_redshift) / \
my_chemistry.a_units
my_chemistry.density_units = mass_hydrogen_cgs # rho = 1.0 is 1.67e-24 g
my_chemistry.length_units = cm_per_mpc # 1 Mpc in cm
my_chemistry.time_units = sec_per_Myr # 1 Myr in s
my_chemistry.velocity_units = my_chemistry.length_units / my_chemistry.time_units
# set initial density and temperature
initial_temperature = 50000. # start the gas at this temperature
# then begin collapse
initial_density = 1.0e-1 * mass_hydrogen_cgs # g / cm^3
# stopping condition
final_density = 1.e12 * mass_hydrogen_cgs
rval = my_chemistry.initialize()
fc = FluidContainer(my_chemistry, 1)
fc["density"][:] = initial_density / my_chemistry.density_units
fc["HI"][:] = 0.76 * fc["density"]
fc["HII"][:] = tiny_number * 0.76 * fc["density"]
fc["HeI"][:] = (1.0 - 0.76) * fc["density"]
fc["HeII"][:] = tiny_number * fc["density"]
fc["HeIII"][:] = tiny_number * fc["density"]
fc["de"][:] = 2e-4 * mass_electron_cgs / mass_hydrogen_cgs * fc["density"]
if my_chemistry.primordial_chemistry > 1:
fc["H2I"][:] = tiny_number * fc["density"]
fc["H2II"][:] = tiny_number * fc["density"]
fc["HM"][:] = tiny_number * fc["density"]
if my_chemistry.primordial_chemistry > 2:
fc["DI"][:] = 2.0 * 3.4e-5 * fc["density"]
fc["DII"][:] = tiny_number * fc["density"]
fc["HDI"][:] = tiny_number * fc["density"]
if my_chemistry.metal_cooling == 1:
fc["metal"][:] = metallicity * fc["density"] * \
my_chemistry.SolarMetalFractionByMass
if my_chemistry.use_dust_density_field:
fc["dust"][:] = metallicity * fc["density"] * \
my_chemistry.local_dust_to_gas_ratio
fc["energy"][:] = initial_temperature / \
fc.chemistry_data.temperature_units
fc["x-velocity"][:] = 0.0
fc["y-velocity"][:] = 0.0
fc["z-velocity"][:] = 0.0
# timestepping safety factor
safety_factor = 0.01
# let the gas cool at constant density from the starting temperature
# down to a lower temperature to get the species fractions in a
# reasonable state.
cooling_temperature = 100.
data0 = evolve_constant_density(
fc, final_temperature=cooling_temperature,
safety_factor=safety_factor)
# evolve density and temperature according to free-fall collapse
data = evolve_freefall(fc, final_density,
safety_factor=safety_factor)
# make a plot of rho/f_H2 vs. T
plots = pyplot.loglog(data["density"], data["temperature"],
color="black", label="T$_{gas}$")
if os.environ.get("METAL_COOLING", 0) == "1":
plots.extend(
pyplot.loglog(data["density"], data["dust_temperature"],
color="black", linestyle="--", label="T$_{dust}$"))
pyplot.xlabel("$\\rho$ [g/cm$^{3}$]")
pyplot.ylabel("T [K]")
pyplot.twinx()
plots.extend(
pyplot.loglog(data["density"], data["H2I"] / data["density"],
color="red", label="f$_{H2}$"))
pyplot.ylabel("H$_{2}$ fraction")
pyplot.legend(plots, [plot.get_label() for plot in plots],
loc="lower right")
if os.environ.get("METAL_COOLING", 0) == "1":
output = "freefall_metal"
else:
output = "freefall"
pyplot.tight_layout()
pyplot.savefig("%s.png" % output)
# save data arrays as a yt dataset
yt.save_as_dataset({}, "%s.h5" % output, data)