-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsave_disk_mass_flow.py
109 lines (89 loc) · 5.71 KB
/
save_disk_mass_flow.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
###############################################
# Save mass and momentum in and around the disk
# for three temperature regimes
###############################################
import yt
from yt import derived_field
import numpy as np
yt.enable_parallelism()
@derived_field(name=("gas","momentum_cylindrical_z"), units="g*cm/s")
def _vertical_flow(field, data):
return data['cell_mass'] * data['velocity_cylindrical_z']
@derived_field(name=("gas","momentum_cylindrical_radius"), units="g*cm/s")
def _radial_flow(field, data):
return data['cell_mass'] * data['velocity_cylindrical_radius']
datasets = yt.load("DD????/DD????", unit_system='galactic')
storage = {}
field_list = ['density','momentum_cylindrical_z','momentum_cylindrical_radius']
for my_storage, ds in datasets.piter(dynamic=False, storage=storage):
dsk = ds.disk([0.5,0.5,0.5], [0,0,1], (20,'kpc'), (1.3, 'kpc'), fields=field_list)
# 10 kpc wider
# buffer1 = ds.disk([0.5,0.5,0.5], [0,0,1], (30,'kpc'), (1.3, 'kpc'), fields=field_list)
# rng = buffer1 - dsk
# Thin ring at edge of disk, 1 radial scale height thick
buffer1 = ds.disk([0.5,0.5,0.5], [0,0,1], (16.5, 'kpc'), (1.3, 'kpc'), fields=field_list)
rng = dsk - buffer1
# 5 kpc taller in each direction
# buffer2 = ds.disk([0.5,0.5,0.5], [0,0,1], (20,'kpc'), (6.3, 'kpc'), fields=field_list)
# slb = buffer2 - dsk
# Thin slabs at top & bottom faces of disk, each 1 vertical scale height thick
buffer2 = ds.disk([0.5,0.5,0.5], [0,0,1], (20, 'kpc'), (0.975, 'kpc'), fields=field_list)
slb = dsk - buffer2
cold_dsk = dsk.cut_region(["obj['temperature'] < 1e4"])
warm_dsk = dsk.cut_region(["(obj['temperature'] > 1e4) & (obj['temperature'] < 1e6)"])
hot_dsk = dsk.cut_region(["obj['temperature'] > 1e6"])
cold_rng = rng.cut_region(["obj['temperature'] < 1e4"])
warm_rng = rng.cut_region(["(obj['temperature'] > 1e4) & (obj['temperature'] < 1e6)"])
hot_rng = rng.cut_region(["obj['temperature'] > 1e6"])
cold_slb = slb.cut_region(["obj['temperature'] < 1e4"])
warm_slb = slb.cut_region(["(obj['temperature'] > 1e4) & (obj['temperature'] < 1e6)"])
hot_slb = slb.cut_region(["obj['temperature'] > 1e6"])
# cold_rad = buffer1.cut_region(["obj['temperature'] < 1e4"])
# warm_rad = buffer1.cut_region(["(obj['temperature'] > 1e4) & (obj['temperature'] < 1e6)"])
# hot_rad = buffer1.cut_region(["obj['temperature'] > 1e6"])
# cold_vrt = buffer2.cut_region(["obj['temperature'] < 1e4"])
# warm_vrt = buffer2.cut_region(["(obj['temperature'] > 1e4) & (obj['temperature'] < 1e6)"])
# hot_vrt = buffer2.cut_region(["obj['temperature'] > 1e6"])
my_storage.result_id = int(ds.basename[-4:])
my_storage.result = [cold_dsk.quantities.total_quantity('cell_mass'),
cold_rng.quantities.total_quantity('cell_mass'),
cold_slb.quantities.total_quantity('cell_mass'),
warm_dsk.quantities.total_quantity('cell_mass'),
warm_rng.quantities.total_quantity('cell_mass'),
warm_slb.quantities.total_quantity('cell_mass'),
hot_dsk.quantities.total_quantity('cell_mass'),
hot_rng.quantities.total_quantity('cell_mass'),
hot_slb.quantities.total_quantity('cell_mass'),
cold_dsk.quantities.total_quantity('momentum_cylindrical_z'),
cold_rng.quantities.total_quantity('momentum_cylindrical_z'),
cold_slb.quantities.total_quantity('momentum_cylindrical_z'),
warm_dsk.quantities.total_quantity('momentum_cylindrical_z'),
warm_rng.quantities.total_quantity('momentum_cylindrical_z'),
warm_slb.quantities.total_quantity('momentum_cylindrical_z'),
hot_dsk.quantities.total_quantity('momentum_cylindrical_z'),
hot_rng.quantities.total_quantity('momentum_cylindrical_z'),
hot_slb.quantities.total_quantity('momentum_cylindrical_z'),
cold_dsk.quantities.total_quantity('momentum_cylindrical_radius'),
cold_rng.quantities.total_quantity('momentum_cylindrical_radius'),
cold_slb.quantities.total_quantity('momentum_cylindrical_radius'),
warm_dsk.quantities.total_quantity('momentum_cylindrical_radius'),
warm_rng.quantities.total_quantity('momentum_cylindrical_radius'),
warm_slb.quantities.total_quantity('momentum_cylindrical_radius'),
hot_dsk.quantities.total_quantity('momentum_cylindrical_radius'),
hot_rng.quantities.total_quantity('momentum_cylindrical_radius'),
hot_slb.quantities.total_quantity('momentum_cylindrical_radius'),]
if yt.is_root():
data_arr = np.zeros((len(storage), 27))
for ds_name, lst in storage.items():
for i in range(data_arr.shape[1]):
data_arr[ds_name, i] = lst[i].v
header = "Disk_Cold_Mgas Ring_Cold_Mgas Slab_Cold_Mgas " + \
"Disk_Warm_Mgas Ring_Warm_Mgas Slab_Warm_Mgas " + \
"Disk_Hot_Mgas Ring_Hot_Mgas Slab_Hot_Mgas " + \
"Disk_Cold_pz Ring_Cold_pz Slab_Cold_pz " + \
"Disk_Warm_pz Ring_Warm_pz Slab_Warm_pz " + \
"Disk_Hot_pz Ring_Hot_pz Slab_Hot_pz " + \
"Disk_Cold_pr Ring_Cold_pr Slab_Cold_pr " + \
"Disk_Warm_pr Ring_Warm_pr Slab_Warm_pr " + \
"Disk_Hot_pr Ring_Hot_pr Slab_Hot_pr"
np.savetxt("disk_mass_flow.txt", data_arr, header=header)