octree total gas masses improperly inversely scaling with nref about yt HOT 11 OPEN

Hi, and welcome to yt! Thanks for opening your first issue. We have an issue template that helps us to gather relevant information to help diagnosing and fixing the issue.

from yt.

Hi! I think this is because the masses shouldn't be smoothed onto the octree, but rather the density, and then that multiplied by the volume. Can you check those results?

That being said, even if that produces the correct answer, this behavior should indeed be flagged in the output so that this doesn't cause issues in the future.

from yt.

Hi! I'm not 100% sure that I did this correctly, so just pasting the volume code I added to the loop in the above script in case I messed up:

    val0 = octree["PartType0","Density"].to('Msun/kpc**3')
    #val = np.sum(octree["PartType0","Masses"].to('Msun').value)
    refined = octree[('index','refined')].astype('bool')
    volume = octree['index','dx'][~refined]*octree['index','dy'][~refined]*octree['index','dz'][~refined]
    val = np.sum(volume.to('kpc**3')*val0)

Assuming that this is right, this is the output I got for the ratios between the total octree mass and the total particle mass:

[1.98063079 2.382016   2.60086731 2.72800857 2.80849894 2.84185174
 2.87547579 2.91535503]

This is obviously not increasing as strongly with nref, but it is, and never seems to be about 1.

from yt.

@DhruvZ this is a bit of a lark but: i'm actually having trouble loading the octree without defining a left edge and a right edge. i wonder if this owes to differing versions of yt or something.

i wonder if loading like this:

octree = ds.octree(ds.domain_left_edge,ds.domain_right_edge,n_ref=64)

makes any difference?

from yt.

I managed to get the same output loading the octree as you described, unfortunately.

from yt.

Based on pixelization_routines.pyx, interpolate_sph_grid_gather (line 1361), I wouldn't expect exact mass conservation from this method: it is calculating the SPH density at cell centers by SPH interpolation to the cell center position, not adding up the mass in each cell. That said, this does not explain why there seems to be a systematic dependence on the octree minimum cell sizes.

Is there any chance that the octree is using the pixelize_sph_kernel_arbitrary_grid function (line 1607) instead? There, I would expect densities to be overestimated in cells that are large compared to particle sizes (see also #4788; I'm working on that one).

from yt.

Hi, thanks for your thoughts! I'm not sure what the most efficient way to check if a was function was called in a package during a script - is there a cleaner way to do it than using pdb and just stepping through?

from yt.

I mean, my idea was the even more basic 'add a print statement to the function'

from yt.

I've been playing around with adding a print statement to the function (and the other related functions), and I'm feel like I'm potentially doing something wrong, but I'm not sure what. I added print statements to the source and reinstalled, but I don't see them out anywhere when running. I do see them where they should be in the compiled code. I also made a new environment for testing this and no change again. This seems to imply that neither of the methods are being called in the script I posted above (either in the original mass-based version or the updated density one).

from yt.

... ok, that's weird. It should be handled by the SPH backend, but maybe the octree stuff has its own routines? I'm not too familiar with those. Do you see the print statements if you add them to the very first thing that gets called? That might tell you if the issue is with the installation or the assumed functions.

I'm trying to trace what gets called here, but that's proving to be a bit tricky. I think

ds.octree(n_ref=nval)

calls the YTOctree class __init__ method (source), which I think should use that pixelization_routines.py function for its SPH smoothing.
On the other hand, the "scatter" option seems to use the CyOctree method for smoothing (source).

from yt.

I just added a print statement to the _determine_fields method in data_objects/data_containers.py and I got many outputs. This certainly suggests that the installation is working correctly. Right now I'm running as in the script I initially posted with gather on, so that is where I also would have expected it to go.

from yt.