_ion_mass and _ion_number_density are not consistent

[claytonstrawn]

Hi all,

I am working on a PR for a new feature allowing a user to specify what metal source they want to use for trident, in case they are comparing multiple codes which do not all retain the same metal species information. However I was struggling to write this PR because I couldn't figure out a clean way to integrate a custom metal function into both priority lists, in _ion_mass and _ion_number_density. But as I was looking into it more, this is partially because having two separate priority lists that look for different fields doesn't really make sense as is.

In my past use of trident, I had written a paper looking at OVI in an ART simulation. I had always only really used the O_p5_number_density field, and never had much need of O_p5_mass. Since ART keeps SNIa and SNII fields separate, I came up with an estimate O abundance that was better than solar, and I added O_nuclei_mass_density fields to the yt ART frontend (yt-project/yt#1981). I didn't realize that this field gets used in the priority list for _ion_number_density but NOT in the nearly-identical list in _ion_mass, and so that field therefore uses solar abundances.

This means, if I take the O_p5_mass field and divide by the volume, it makes something very different from O_p5_density, even though as a user I would naively assume them to be the same. The difference between them can be pretty dramatic:

import yt,trident,agora_analysis

snap = agora_analysis.AgoraSnapshot('art',4.0)
snap.load_snapshot()
trident.add_ion_fields(snap.ds,['O VI'])

def ion_density_from_mass(field,data):
    if isinstance(field.name, tuple):
        ftype = field.name[0]
        field_name = field.name[1]
    else:
        ftype = "gas"
        field_name = field.name
    ion = field_name.split('_density')[0]
    volume = data['gas','cell_volume']
    mass = data['gas','%s_mass'%ion]
    return mass/volume
snap.ds.add_field(('gas','O_p5_density_from_mass'), function=ion_density_from_mass, units='g/cm**3',
                    sampling_type='cell')

p1 = yt.ProjectionPlot(snap.ds,0,('gas','O_p5_density'),center = snap.center,width = 2*snap.Rvir)
p1.set_zlim(('gas','O_p5_density'),1e-10,1e-5)
p1.show()
p2 = yt.ProjectionPlot(snap.ds,0,('gas','O_p5_density_from_mass'),center = snap.center,width = 2*snap.Rvir)
p2.set_zlim(('gas','O_p5_density_from_mass'),1e-10,1e-5)
p2.show()

You can ignore the "agora_analysis" stuff, that's just basically wrappers for yt.load and containing the zoom-in center info, etc. This same effect will be true for any ART dataset, because my changes are in the frontend. But really, this will happen to any user who defines O_nuclei_mass_density without defining O_nuclei_mass or vice versa. Which I imagine could be pretty common.

I guess I just want to ask if there's any interest in updating the logic of ion_balance to have only one priority list, and have ion_mass be derived from ion_number_density, like ion_density currently is. This would force these fields to be consistent.

I'm happy to write the PR to do that if so, but wanted to raise this as an issue regardless.

[chummels]

Hi Clayton. Thanks for the writeup. I can confirm that this appears to be an issue. I recall that when we wrote this originally, all the ion fields (mass, density, abundance, and number density) were directly tied to each other, such that this wasn't a problem. I think they potentially got tweaked with your PR from 4 years ago: PR #81 , which I think is the last to touch this functionality. I'll go through things and try to figure out the desired behavior (ie which field is correct). When I test this locally on the Mg II field, there is a fixed multiplicative offset between the two fields (density and density_from_mass) of 0.0433, which may be a hint, but I need to dig into a bit to figure out where things are going wrong.

[claytonstrawn]

Wow! I had completely forgotten my part in writing this for #81 back in 2019! Well, in looking at some of @brittonsmith's comments there, I think the main reason I added this second sequence of priority checks is there is because we were not confident that SPH codes would have a "volume" we could multiply by density to get mass. One thing I've been doing with SPH codes in AGORA is letting yt's internal SPH handling do the heavy lifting, as each gas element surely has a mass and density.

effective_volume = data['gas','mass'] / data['gas','density'] 
ion_mass = data['gas','O_p5_density'] * effective_volume

I don't think it necessarily matters whether effective_volume really represents a physical volume or just a sort of inverse SPH smoothing, either way it just says, "whatever the relationship is for a gas parcel's density to its mass, use that same relationship between the ion's density and mass". Would this be a good enough solution?