[WIP] [bugfix] yt 4.0 deposit bug fix

[AshKelly]

PR Summary

I've added some support such that covering_grid and arbitrary_grid in the demeshed branch can accept,

arb[('deposit', 'PartType0_mass')]

as well as,

arb[('gas', 'mass')] 

fixes #2003

This isn't currently working as expected

PR Checklist

• [x] Code passes flake8 checker
• [ ] NA New features are documented, with docstrings and narrative docs
• [ ] NA Adds a test for any bugs fixed. Adds tests for new features.

[ngoldbaum]

Looks like some of the arbitrary_grid tests are failing.

Can you also add tests for this use case? I'm pretty sure the changes you've made will work for covering_grid and aribtrary_grid out of the box, but best to explicitly test that.

Thanks for sending in the patch and working on this!

[AshKelly]

hm, so this doesn't work. We get to a point where if you pass in ('PartType0', 'mass') or ('deposit', 'PartType0_mass') that we can't tell if we should read from the file or SPH interpolate.

The split fields decomposes saying we need, [('PartType0', 'Coordinates'), ('PartType0', 'Masses')] and then we don't know whether to read the Coordinates and Masses or try to interpolate as field[0] is the SPH particle type. Obviously we should ever interpolate coordinates, but in some instances people may want the Masses to be interpolated.

[ngoldbaum]

When you get passed ('deposit', 'PartType0_mass') you could short-circuit directly to the deposition field machinery, and get the underlying particle field data from the data object the covering grid uses for data access (grid._data_source).

[AshKelly]

It works (I think) for all other particles except PartType0 but I don't think I currently understand the flow of this section of code to figure this out after all. I'll have to come back to it later and try to understand more what is happening with the stars so I can replicate it for the gas

[ngoldbaum]

OK. If you need a hand just let me know, I'm pretty sure I originally wrote the relevant code.

[AshKelly]

yeah, so, what seems to happen is the ('deposit', 'PartType0_mass') comes in and then this goes into generate_fields which then sends ('PartType0', 'mass') and ('PartType0', 'Coordinates') through. I can't seem to get a handle on how to distinguish against that and when somebody directly asks for ('PartType0', 'mass')

[AshKelly]

so this is my test script

import yt

fname = 'GadgetDiskGalaxy/snapshot_200.hdf5'
ds = yt.load(fname)

center = [31995.63476562, 31473.6640625, 28969.88671875]
arb_center = ds.arr(center, 'code_length')
left_edge = arb_center-ds.quan(100, 'kpc')
right_edge = arb_center+ds.quan(100, 'kpc')

arb = ds.arbitrary_grid(left_edge, right_edge, dims=[128, 128, 128],
                        field_parameters={'center': center})
ds.index

print("#### Star mass deposit")
arb[('deposit', 'PartType4_mass')].in_units('Msun')

print("#### Gas mass deposit")
arb[('deposit', 'PartType0_mass')].in_units('Msun')

Essentially the call from fields/particle_fields.py which I have added a print for makes a call which to me is indistinguishable from someone asking for ['gas', 'mass']

[matthewturk]

@AshKelly is this still something we should address?

[AshKelly]

yeah, so I'm still stuck with this. In summary, my issue is that we want to intercept ("gas", "XXX") and as its an SPH field we smooth it, but then it was suggested that ("PartType0", "XXX") should also get SPH smoothed.

The issue is when ('deposit', 'PartType0_mass') is passed in, it gets decomposed into ('PartType0', 'mass') and ('PartType0', 'Coordinates') which then revisit the initial function. So then I try to SPH smooth the ('PartType0', 'mass') which is wrong.

An easy-ish fix is to just SPH smooth ("gas", "XXX") I guess or to make SPH smoothing a thing e.g, ('smooth_deposit', 'PartType0_mass')

All opinions welcome on this :)

[matthewturk]

This might be a place where we need to break the similarity between mesh and particle codes, and make it clear that we are explicitly smoothing something. If we agree to that -- and it might need to come up on yt-dev -- then I think the final option you propose would be good.

[munkm]

Hey @AshKelly do we want to add the 4.0 label for this PR?

[AshKelly]

Yeah, this a big thing we should fix. This is related to the whole gas field name API stuff.

[AshKelly]

TLDR

We wired up the covering_grid and arbitrary_grid so that if they detect the particle type PartType0 in the get_data function then we short-circuit and go into some SPH deposition.

[see, L686 of yt/data_objects/construction_data_containers.py in this branch]

        if(hasattr(self.ds, '_sph_ptype') and len(part)>0):
            is_sph_field = self.ds._sph_ptype in part[0]

but, sometimes people want to use the older machinery of depositing gas and pretending it is not an SPH field, e.g. if we care about total mass, arb[('deposit', 'PartType0_mass')]. Unfortunately, the method get_data recursively identifies that it needs the fields ('PartType0', 'Coordinates') and ('PartType0', 'mass') and calls itself recursively with those fields.

Then we go "ah! this is an SPH field, let's SPH smooth!" because we have forgotten the original function call. This means we end up trying to SPH smooth, but the user didn't ask for that.

So, I'm happy to work on this, but does anyone have any clever suggestions how I may fix this? If not, I'll just go along on my own and see if I can dig something out which works.

Pinging @munkm @neutrinoceros @matthewturk (for both suggestions and to make sure we keep this alive and fix before release)

[stonnes]

Hi all, I am making my way through our PRs, and noticed that things were moving along here but seem to have hit a snag on SPH smoothing of particles. @AshKelly were you ever able to figure that out? @munkm @neutrinoceros @matthewturk you were all pinged to see if you had ideas on the problem detailed above, this is just an extra nudge to make sure it doesn't fall off your radar! Thanks!

[AshKelly]

I think this is certainly something we should fix, however it's stagnant because I'm lost on how to fix it. I'm happy to zoom with someone and sort of discuss, brainstorm and hopefully fix, but I don't understand the field system well enough to fix this alone.

[brittonsmith]

For posterity, I believe this is related to PR #2986.

[munkm]

@AshKelly I definitely think a zoom bugfix session is a good idea. Maybe we can poll people on slack and the mailing list to see who wants to join. I am happy to join but maybe another person (@jzuhone @Xarthisius or @matthewturk maybe?) would be interested too?

[matthewturk]

@AshKelly I've thought through this a bit, and I have an opinion, but I'm not sure how strongly I hold it.

I think that if we're asking for density, specifically, it should put the density down -- and that means to smooth it. I think that asking for a particle field name, for instance Density, would be fine to add as a non-smoothed field. But I don't completely know the right way to implement this, unless we also did something like either check if it's an alias field (thus smoothed) or to add, say, a Particle Union that maps the PartType0 to some other field name. That would mean that if you accessed, say, ("PartType0", "density") you'd get smoothed, but if you made a new alias to "PartType0_unsmooth" (or something) then when you accessed ("PartType0_unsmooth", "density") it would just jam the values in.

I think that density is maybe the bad field example (that I chose!) because stuff like "Temperature" might be the ones you wouldn't want smoothed (plus "Density" is also calculable by doing mass/volume, and volume here will be a free paramter.)

Anyway, kinda what I'm thinking is that we probably don't want to support raw deposition of particles as opposed to smoothed, because doing so can be ill-defined, and requires more attention to the intensive/extensive properties of the field. In the cases where you do want to do that, it might be better suited to a raw binning via something like Profile3D, although we don't have any docs that say that.

[AshKelly]

@matthewturk I'm happy with both suggestions. Dropping the support, but also allowing some kind of override name ("PartType0_unsmooth", "density").

I was actually trying to implement the latter but due to the recursive nature of the get_data function the initial call gets washed out. If someone asks for ("PartType0_unsmooth", "density") then get_data says "Ok, I need positions and densities for PartType0". It then recursively calls itself with one of the calls asking for ("PartType0", "density"). I can't distinguish this call from a user asking for SPH smoothing of ("PartType0", "density") so the code returns something of the wrong shape and then everything breaks.