[AstroPy] Implementing a versatile natural unit system with astropy?

Wed Jan 24 19:17:40 EST 2024

Hi Nick,

I was hoping for you to chime in, since I've thought of submitting this as a feature request on plasmapy. But I haven't been keeping up with plasmapy updates due to a few recent side quests and am not sure how the simulation module wants to progress over there.

What you said was one of my personal code iterations too - to use objects like [conversion factors in this documentation](https://github.com/lanl/vpic/blob/master/doc/vpicUnits.pdf) for the Vector PIC code. I went as far as writing a Quantity subclass (SimQuantity) with normalization attributes and methods likeSimQuantity.code_to_user (returning a dimensional quantity) and SimQuantity.user_to_code (returning a dimensionless quantity).

But I realized this level of abstraction may not be necessary when we can define the conversion factors with astropy.units.def_unit like

L_factor = u.def_unit("L_factor", 10 * u.km)
T_factor = u.def_unit("T_factor", 0.5 * u.ns)
Q_factor = u.def_unit("Q_factor", c.si.e)
M_factor = u.def_unit("M_factor", c.si.m_p)

... and populate these factors with a method like astropy.units.natural.set_system(system) where system is one of those in Section 3 of the linked document, which does

u.add_enabled_units([L_factor, T_factor, Q_factor, M_factor])

... so that they are searchable via UnitBase.find_equivalent_units(). Now, everyone (data scientists and simulators) can all just use u.Quantity for their own purposes. For example,

Data scientist:
>>> obs_data = load_obs_data() # magnetic field, electric field, ion cyclotron frequency, etc
>>> q = some_data_science_calculations(obs_data)
... then, to know what q is in ion scales>>> u.natural.set_system("ion_scale", ...) # inject parts of obs_data here>>> q.to("L_factor")
<Quantity 0.02 ion_scale L_factor>

Simulators:
>>> u.natural.set_system("ion_scale", ...) # inject motivating obs_data here
>>> p = 0.5 * u.natural.L_factor
>>> q = run_simulation(p)
>>> q.to("km") # to see what q is in si
<Quantity 10. km>

In these examples, no one has to worry about keeping track of conversion factors themselves, while working on their units of choice and only having intuitive objects storing the dimensional data. And we don't have to subclass astropy.units.Quantity at all. I suppose plasmapy could host such a module like astropy.units.natural (i.e., plasmapy.natural_units; note that it doesn't have to be under plasmapy.simulation). But this will result in loss of generality, because all of the above are not specific to plasma physics.

-- T

On Wednesday, January 24th, 2024 at 10:51 PM, Nick Murphy via AstroPy <astropy at python.org> wrote:

> Thank you for bringing this up, and also the very insightful responses!
>
> Over in PlasmaPy, I've started to work on classes for representing normalizations of different systems of equations, building on astropy.units. The idea is to create normalization objects which contain attributes with the normalizations for length, time, and other physical types. My goal is to enable operations like:
>
> dimensionless_object = dimensional_object / normalizations_object
> dimensional_object = dimensionless_object * normalizations_object
>
> where dimensionless_object and dimensional_object could be single or compound Quantity objects. One of the main purposes would be to facilitate comparisons between simulations and physical systems.
>
> I've been meaning to get back to an old [pull request to add an MHDNormalizations prototype](https://github.com/PlasmaPy/PlasmaPy/pull/1029) as a first step, based on the normalizations for the MHD equations. The idea is to provide MHDNormalizations with three quantities representing different types (e.g., length, magnetic field, & number density), and then (for example) get the normalization for time from the time attribute. This is still a work in progress, so I'd be happy to talk more about this and learn more from you about your perspective.
>
> Thank you again!
>
> On Wed, Jan 24, 2024 at 2:38 PM John ZuHone <jzuhone at gmail.com> wrote:
>
>>> On Jan 24, 2024, at 2:19 PM, Tien via AstroPy <astropy at python.org> wrote:
>>>
>>> From a user perspective, if I only care about doing simulations, then unyt seems like a good library. But then if I want to use astropy and its other utilities, I'll have to manually handle (kind of redundant) conversions between unyt types and astropy types. Wouldn't it still be good to have native support in astropy?
>>
>> Yep, I agree that it would be a nice thing to have. I use both unyt and astropy.units in my work, usually in different contexts (and also where the former isn’t guaranteed to be installed), but they are both excellent implementations. If there’s any way I could help with an astropy implementation I’d be happy to.
>>
>> The main reason the connecting tissue between the two is somewhat kludgey is simply because they use different machinery under the hood—unyt uses sympy to handle all of its symbolic unit operations whereas astropy has rolled its own implementation, so we wouldn’t be able to marry them trivially.
>> _______________________________________________
>> AstroPy mailing list
>> AstroPy at python.org
>> https://mail.python.org/mailman/listinfo/astropy
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/astropy/attachments/20240125/258775b7/attachment-0001.html>