calculate the mass flux in yt
Dear all, I have a question on calculating the mass flux across surfaces in yt, particularly why the flux values differ between two methods. The surface, for example, is a circle a few kpcs above the central XY plane, and I'm calculating the mass flux in the perpendicular (velocity_z) direction. A simplified code snippet looks like this. ds = yt.load('Enzo_AMR_output') box_length_kpc = ds.domain_width[0].in_units('kpc').value # constant, simulation box length in kpc *Method 1.* Create a 2D yt.slice object and calculate the flux manually (faster). slice = ds.r[:, :, 0.5 + 2/box_length_kpc] # selecting a slice at z=2 kpc (relative to domain center) # creating a boolean mask: a circle with a radius of 10 kpc radius_mask = np.sqrt((slice['gas','x'].in_units("kpc").value-box_length_kpc/2.)**2 +\ (slice_plus['gas','y'].in_units("kpc").value-box_length_kpc/2.)**2) <= 10 # mdot = rho * v dot dA, where dA = dx * dy mass_flux = np.sum((slice['gas', 'density'][radius_mask_] *\ slice['gas', 'velocity_z'][radius_mask] *\ slice['gas', 'dx'][radius_mask] *\ slice['gas', 'dy'][radius_mask])).in_units("Msun/yr") *Method 2.* Create a yt.surface object and use the calculate_flux() function. # defining a disk object of radius 10 kpc, then selecting the upper surface where z=2 kpc (relative to domain center), and calculate the mass flux my_disk = ds.disk(ds.domain_center, [0.0, 0.0, 1.0], 10 * kpc, 2.02 * kpc) disk_upper_surf = ds.surface(my_disk, ("gas", "cylindrical_z"), 2) mass_flux = disk_upper_surf.calculate_flux(("gas", "velocity_x"),("gas", "velocity_y"),("gas", "velocity_z"),("gas", "density")).in_units("Msun/yr") The two methods always return different results (off by ~48% for this selected input), but I wonder why. Another sanity check I did was to print out the cell numbers selected in each method. Method 2 (marching cubes, triangular vertices) consistently selected 2 times the cell numbers in method 1 (slicing). Could anyone point out what might be wrong with the way I implemented either method? Thanks so much! Best, Jingyao -- Jingyao Zhu (she/her/hers) Columbia University | Dept of Astronomy
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
Hi Nina, Thanks so much for the suggestions and the helpful tests! For method 2, I defined the disk to have a slightly larger height (e.g., 2.02 kpc) for selecting an isocontour surface of cylindrical_z = 2.0 kpc to avoid having all zeros. If I printed the selected values, it would be an array of [2.0, 2.0, ..., 2.0], because the ds.surface sampling method (marching cubes) interpolated to that exact planar value. So methods 1 and 2 will not have the same data selection by default (one is grid-based, the other from interpolation). I still wonder why these two surface-based methods would return different fluxes. Thanks for suggesting the volume-based solution using ds.include_above/below, I will try this! For the yt developers/community users, I would still love to hear your insights on the surface.calculate_flux method: It seems very useful for, e.g., calculating fluxes through curved surfaces like the side of a disk or a sphere, but since this test on a flat surface doesn't pass, I'm wondering how to interpret the results. Thanks again! Jingyao On Wed, Apr 26, 2023 at 8:08 AM Nina Akerman <akermannia@gmail.com> wrote:
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
-- Jingyao Zhu (she/her/hers) Columbia University | Dept of Astronomy
participants (2)
-
Jingyao Zhu -
Nina Akerman