[yt-dev] Metallicity units

[Nathan Goldbaum]

Hi all,

I've updated the metallicity units PR:
https://bitbucket.org/MatthewTurk/yt/pull-request/38/fixing-cgs-unit-conversions-for-angle-and/diff

I ended up deciding that metallicity doesn't really make much sense as
a dimension in the unit system.  Instead, I've made it so the
`code_metallicity` and `Zsun` units are dimensionless.

Here's an example run using the tip of the PR which illustrates why
this is a nice choice:

First, load the data:

>>> pf = load('HiResIsolatedGalaxy/DD0044/DD0044')
>>> dd = pf.h.all_data()

Get the mean metal mass for all particles in the simulation:

>>> mean_metal_mass = (dd['metallicity_fraction']*dd['particle_mass']).mean()
>>> print mean_metal_mass
4.13009892602e-09 code_mass*code_metallicity
>>> print mean_metal_mass.in_cgs()
8.21517977375e+33 g

If I had run this using a version with a metallicity dimension, I
would have had to coerce the units to be grams, even though code
metallicity is really just metallicity fraction. I think getting a
metal mass or density like this will be an extremely common operation,
and it's something that a metallicity unit needs to handle naturally.

Here's another example where I get the mean metallicity of all
particles in the simulation:

>>> mean_metallicity = dd['metallicity_fraction'].mean()
>>> mean_metallicity.in_cgs()
0.00503378562375 dimensionless
>>> mean_metallicity
0.00503378562375 code_metallicity
>>> mean_metallicity.in_units('Zsun')
0.246633298567 Zsun

In principle code_metallicity can have a non-unit conversion to a
dimensionless fraction, but enzo writes out particle metallicities as
mass fractions, so that is unnecessary here.

As another concrete example why this way is better, here is the way I
needed to define the metal_density field for a Tipsy SPH dataset:

@derived_field(name = "metal_density", units = "g/cm**3")
def metal_density(field, data):
    ret = data['metallicity_fraction']*data['density']
    ret.convert_to_units('g/cm**3*code_metallicity')
    ret = np.array(ret)
    return data.pf.arr(ret, input_units='g/cm**3')

If metallicity weren't a base dimension, I wouldn't need to coerce the
units to g/cm**3, the conversion to CGS would handle that
automatically.

Hope that explains my position clearly.  I'd love to hear comments,
concerns, or alternate approaches - this is definitely tricky.

-Nathan