Hi Jingyao,
Regarding the method 2, I have never worked with ds.surface, but looking at yt wiki I’m not sure it is what you need here. As a sanity check, I asked to print cylindrical_z of disk_upper_surf cells, and I think it’s all zeros, since I never got printed anything:
for i in range(len(disk_upper_surf['cylindrical_z'])):
if disk_upper_surf['cylindrical_z'][i] < 0:
print('gotcha’)
for i in range(len(disk_upper_surf['cylindrical_z'])):
if disk_upper_surf['cylindrical_z'][i] > 0:
print('gotcha')
So with ds.surface you select exactly the planar cells.
The way I usually select data is with include_above/include_below:
my_disk = ds.disk(ds.domain_center, [0.0, 0.0, 1.0], 10 * kpc, 2.02 * kpc)
disk_upper_surf = my_disk.include_above('cylindrical_z', 0, 'kpc')
Could you please check if it works? Then you I think it would be good to compare the results method 1 and method 2 give on exactly the same data selection, so mass_flux = np.sum(disk_upper_surf[‘density'], ...)
Finally, I haven’t used calculate_flux either, frankly, I didn’t know it even existed. Again, from the wiki it’s not quite clear how it works, especially if you’re interested in the flux specifically from velocity_z, not from all the velocities. Here is my (probably overly complicated) solution to this problem that I often use in my work:
from yt.units import kpc #(we’ll need it)
def _momentum(field, data):
return data['cell_mass'] * data['velocity_cylindrical_z’]
ds.add_field(('gas', 'momentum'), units='kpc * Msun/yr', function=_momentum, sampling_type='cell’)
my_disk = ds.disk(ds.domain_center, [0.0, 0.0, 1.0], 10 * kpc, 2.02 * kpc)
disk_upper_surf = my_disk.include_above('cylindrical_z', 0, 'kpc')
inflow = disk_upper_surf.include_below('momentum_z', 0)
outflow = disk_upper_surf.include_above('momentum_z', 0)
mass_inf = inflow.sum('momentum_z') / (height * kpc)
mass_outf = outflow.sum('momentum_z') / (height * kpc)
The idea here here is that mass flux is sum(m * v_z) / z
So I first create a new cell type that is momentum p = m * v. inflow and outflow allow me to get negative and positive mass flux values separately, which I do with mass_outf and mass_outf, where height is the vertical height, which in your case is 2.02 kpc (from 0 to the height in your my_disk definition). Btw, why is it 2.02 not just 2?
Mass flux is then simply mass_inf + mass_outf
Hope it helps!
Let me know,
Nina