Some updates and additions to astrophysics-related fields #2245

jzuhone · 2019-04-17T15:01:29Z

PR Summary

This PR updates a number of the astronomy/astrophysics-specific fields in yt. In detail:

Wherever possible, change "cell_mass" to "mass" so that fields are not only defined for grid-based datasets. Similarly, if possible, do not use "cell_volume" for fields which are meant to be general.
Add a computation for the number density fields for free electrons, protons, and helium nuclei for when there are no species fields, assuming primordial abundances and full ionization.
Redefine "number_density" and "mean_molecular_weight" fields so that they can be calculated even if no species are defined, by assuming a constant mean molecular weight and primordial abundances.
Add an "optical_depth" field which depends on the "El_number_density" field. Redefine the two SZ fields in terms of this field (otherwise they just assume full ionization, which may not be true).
Add a "velocity_los" field for projections, which computes the velocity component along a line of sight, on-axis or off-axis. Use it in the "sz_kinetic" field.
Add a "magnetic_field_los" field for projections, which computes the magnetic field strength component along a line of sight, on-axis or off-axis. This is used for the Faraday rotation measure field, which is also added here.
Add an alias for the "particle_gpot" field in the FLASH frontend.
Move the "entropy" field to astro_fields, since its field definition is very astro-specific (and in fact it's very high-energy-astro-specific). Make it depend on the "El_number_density" field (otherwise it just assumes full ionization, which may not be true).
Remove chandra_emissivity field. It's not well-documented and it is not clear at all what it represents (i.e. what is the energy band? what is the instrument assumed?), so it's misleading.
Remove extra whitespace and parentheses in various places
Set up a function which computes a default mean molecular weight for frontends that may use one for defining things like temperature and number density but don't have one. This default mu is now used in the athena, athena_pp, and gadget frontends unless the user passes in a different value when loading the dataset. NOTE: This default mu assumes full ionization.
Made sure that the flash, enzo, and gamer frontends get to define fields like "number_density" the way they want.
Fixed a bug which prevented magnetic fields from being defined properly in SPH frontends.
The "cell_mass" and "courant_time_step" fields should have sampling_type="cell".
Add species fields to the Arepo frontend.
Fixed some bugs created in off-axis projections

PR Checklist

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

jzuhone · 2019-07-05T19:21:07Z

After rebasing this PR I'd like to come back to it. The major unresolved issue is what I said above (see the full text above for context):

There is one subtle inconsistency in this setup at the moment. The field definitions I wrote for the "number_density' and "mean_molecular_weight" fields assume primordial abundances and zero ionization. This is what we do in a number of the SPH frontends when we need to compute the temperature but there are no species fields to compute the molecular weight. But defining an free electron number density field really only makes sense if at least the gas is partially ionized, so in that case I've assumed full ionization. But there's obviously an inconsistency there. It would be easy to just define the number density fields here and the temperature fields in the SPH frontends to assume full ionization instead if there are no species defined, but this is potentially a significant backwards-incompatible change for the latter.

The possible resolutions I can see are:

Have the default for both species fields and the mean molecular weight be zero ionization
Have the default for both species fields and the mean molecular weight be full ionization
Have this set up in a frontend-dependent way, perhaps by using parameters

brittonsmith

@jzuhone, these code improvements all look good to me. My only suggestion is to perhaps create an add_los_field function in vector_operations.py to generate line of sight fields from vector fields.

On the mean molecular weight issue, my inclination is toward assuming full ionization. As you mention, the neutral assumption leaves you with zero electron density, requiring the extra inconsistency to model that simply. This might be just me, but my impression is that simulations that focus on neutral gas are more likely to have the additional fields to compute mmw correctly, whereas simulations of higher temperature phenomena are more likely to follow just density.

In Grackle, we also include the effects of metals by assuming a mean molecular weight of 16 and using the metal density field. That number is about what you get for assuming fully ionized metals with a solar abundance pattern.

jzuhone · 2019-07-08T13:35:50Z

@brittonsmith I like that idea regarding vector operations.

Anyone who works with SPH data have any concerns about changing the zero ionization assumption? @qobilidop or @AshKelly? Note that this only affects data which doesn't have any abundances defined.

Relevant code is here, particularly the last 10 lines or so:

yt/yt/frontends/gadget/fields.py

Lines 70 to 92 in 33dcf44

    
           def setup_gas_particle_fields(self, ptype): 
        
               if (ptype, "ElectronAbundance") in self.ds.field_list: 
        
                   def _temperature(field, data): 
        
                       # Assume cosmic abundances 
        
                       x_H = 0.76 
        
                       gamma = 5.0/3.0 
        
                       a_e = data[ptype, 'ElectronAbundance'] 
        
                       mu = 4.0 / (3.0 * x_H + 1.0 + 4.0 * x_H * a_e) 
        
                       ret = data[ptype, "InternalEnergy"]*(gamma-1)*mu*mp/kb 
        
                       return ret.in_units(self.ds.unit_system["temperature"]) 
        
               else: 
        
                   def _temperature(field, data): 
        
                       # Assume cosmic abundances 
        
                       x_H = 0.76 
        
                       gamma = 5.0/3.0 
        
                       if data.has_field_parameter("mean_molecular_weight"): 
        
                           mu = data.get_field_parameter("mean_molecular_weight") 
        
                       else: 
        
                           # Assume zero ionization 
        
                           mu = 4.0 / (3.0 * x_H + 1.0) 
        
                       ret = data[ptype, "InternalEnergy"]*(gamma-1)*mu*mp/kb 
        
                       return ret.in_units(self.ds.unit_system["temperature"])

jzuhone · 2019-07-08T22:15:12Z

Would also love any comment from @chummels

brittonsmith · 2019-07-09T08:25:43Z

@jzuhone, I think Cameron is busily preparing to leave on a long trip at the end of the week. It's likely he's not going to have time to comment here.

AshKelly · 2019-07-09T09:35:37Z

yeah, I think zero ionization is appropriate if there are no other fields, e.g electron fraction. Hopefully, if people don't explicitly store electron fraction or temperature then they aren't too worried about it

jzuhone · 2019-07-09T12:59:48Z

@AshKelly the proposal here currently is to assume full ionization unless specified.

Unfortunately it’s not true that if people don’t store those things that they don’t care (although in theory it should be). For most applications in X-ray astrophysics, full ionization is assumed and so you need a mean molecular weight which reflects that to get electron number density and temperature, both of which are needed for emission and other calculations.

So we need to decide if we want to assume full ionization by default (a decidedly higher-temperature astrophysics bias) or zero ionization (a lower temperature general physics bias). Or if we want to have some kind of switch that does either.

matthewturk · 2019-07-12T00:25:23Z

@jzuhone This is a really interesting question. I'm not sure I have the right answer, but I would propose that since we're talking about a default, we make room for different codes and frontends to set it, and we also identify how to set it -- probably in something less fragile than field_parameters.

I think that if we can't deduce it, then it's safe for this field to assume full ionization. Does that make sense?

jzuhone · 2019-07-12T13:47:01Z

Ok, I made a flurry of changes which I think settles the question regarding how to set mean molecular weights. Please have a look at the latest changes starting with 3e412fd and leave some feedback.

brittonsmith · 2019-07-16T15:48:14Z

yt/fields/vector_operations.py

+                       for ax in axis_order]
+
+    def _los_field(field, data):
+        axis = data.get_field_parameter("axis")


Quick question, does axis need to be checked to see if it's a unit vector?

It could be either a unit vector or a single integer for an axis-aligned l.o.s. Do you have a suggestion for a check?

brittonsmith · 2019-07-16T15:52:17Z

yt/frontends/athena_pp/fields.py

-            else:
-                mu = 0.6
-            return (data["gas","pressure"]/data["gas","density"])*mu*mh/kboltz
+            return (data["gas","pressure"]/data["gas","density"])*data.ds.mu*mh/kboltz


Looks like above you're setting self.parameters["mu"]. Do you need to use that here?

Ah, I didn't do this right. Thanks for the find again.

brittonsmith · 2019-07-16T15:53:49Z

yt/frontends/enzo/fields.py

+                ("gas", "mean_molecular_weight"),
+                sampling_type="cell",
+                function=_mean_molecular_weight,
+                units=unit_system["number_density"])


Are these the correct units?

No, you are right, I'll fix that. Thanks for the catch.

brittonsmith · 2019-07-16T15:57:08Z

yt/frontends/flash/fields.py

+                           function=_number_density,
+                           units=unit_system["number_density"])
+
+        if "ye" in self.field_list:


I'm not that familiar with the FLASH format. Can you get both this condition and the one above (if ("flash", "abar") in self)? If so, you could have one number_density field overwriting another. Should this be an elif?

brittonsmith · 2019-07-16T16:06:24Z

yt/utilities/chemical_formulas.py

+         (2.0 * _primordial_mass_fraction["H"])
+    mu += ChemicalFormula("He").weight / \
+          (3.0 * _primordial_mass_fraction["He"])
+


This never returns mu.

yt/utilities/chemical_formulas.py

yt/fields/species_fields.py

brittonsmith

I left one last minor comment, but this looks good to me. Cameron is away on a hike for the next 3+ months and I imagine Nathan is pretty busy with other things. @matthewturk, do you think you could take a quick look at this? I think two approves would be sufficient.

Also, John, we have a bot that prevents WIP PRs from being merged, so please go ahead and remove that from the title when you're ready to get this in.

brittonsmith · 2019-07-23T12:06:53Z

yt/frontends/flash/fields.py

+                Na_code = data.ds.quan(Na, '1/code_mass')
+                return data["flash", "dens"]*Na_code/data["gas", "mean_molecular_weight"]
+
+            self.add_field(("gas", "number_density"),


Minor, but you could indent this self.add_field call one to the left and remove the identical call inside the if statement above.

jzuhone · 2019-08-08T14:20:44Z

Since this touches a lot of stuff, we should probably enforce three approvals here.
@qobilidop, @AshKelly, @munkm, anyone?

jzuhone · 2019-08-14T17:06:03Z

Anyone want to take a peek at this?

matthewturk

Impressive and detailed! Good stuff.

doc/source/analyzing/fields.rst

yt/data_objects/construction_data_containers.py

yt/fields/astro_fields.py

matthewturk · 2019-08-16T16:16:45Z

I see no reasons to hold off on merging this.

jzuhone · 2019-08-16T16:22:44Z

@matthewturk fine by me! Thanks!

munkm

This looks really great! I really like the line of sight field feature. I have a few minor things I've commented on, mostly docs clarifications.

doc/source/analyzing/fields.rst

munkm · 2019-08-17T00:55:03Z

doc/source/analyzing/fields.rst

+
+At this time, this functionality is enabled for the velocity and magnetic vector
+fields, ``("gas","velocity_los")`` and ``("gas","magnetic_field_los")``. The 
+following fields built into yt make use of these line-of-sight fields:


If a user wants to have a line of sight field for a different vector field, is there a way to do that?

Not yet, but in fact, none of the vector field operations (like magnitude, etc) are available for any other fields than these two. We should probably expose this in a separate PR somehow.

doc/source/examining/loading_data.rst

yt/fields/xray_emission_fields.py

jzuhone · 2019-08-19T14:46:31Z

@munkm I think I've addressed everything, but let me know if there is anything else. Thanks for the comments!

munkm

Thanks for the updates @jzuhone! The clarifications you added were really helpful!

jzuhone requested review from ngoldbaum, brittonsmith and chummels April 17, 2019 15:01

jzuhone force-pushed the astro_fields branch 4 times, most recently from 1dc4392 to 6c811d7 Compare April 24, 2019 12:54

jzuhone force-pushed the astro_fields branch from 6c811d7 to ec9d051 Compare May 29, 2019 20:23

jzuhone force-pushed the astro_fields branch 2 times, most recently from 5fe401d to da34055 Compare July 5, 2019 19:01

brittonsmith approved these changes Jul 8, 2019

View reviewed changes

jzuhone force-pushed the astro_fields branch from 3269c3a to 82018d5 Compare July 14, 2019 12:50

brittonsmith reviewed Jul 16, 2019

View reviewed changes

yt/utilities/chemical_formulas.py Show resolved Hide resolved

yt/fields/species_fields.py Show resolved Hide resolved

jzuhone force-pushed the astro_fields branch 2 times, most recently from a13760c to e7c14a4 Compare July 22, 2019 19:00

brittonsmith reviewed Jul 23, 2019

View reviewed changes

jzuhone added 9 commits August 6, 2019 14:25

Only add this field if we have to

4619831

This is the correct way to make sure all of these fields are set up

55970e5

Add these fields for Arepo testing

eca2731

Avoid conflict

c6a7156

This goes in here instead

4c1d6e8

This is the correct logic here

1791d8a

Make sure we can do this for on-axis projections

39fb9a8

Final fixes for filtered particles and off-axis projections

2c5f265

A note of warning

9ccbf5a

jzuhone force-pushed the astro_fields branch from ec6651c to 9ccbf5a Compare August 6, 2019 18:27

jzuhone added 3 commits August 6, 2019 14:32

Attempt to get appveyor tests to pass

8eef30e

Bugfix

4235521

Bugfix

11e3d09

brittonsmith approved these changes Aug 8, 2019

View reviewed changes

matthewturk approved these changes Aug 15, 2019

View reviewed changes

doc/source/analyzing/fields.rst Show resolved Hide resolved

doc/source/analyzing/fields.rst Show resolved Hide resolved

yt/data_objects/construction_data_containers.py Outdated Show resolved Hide resolved

yt/fields/astro_fields.py Show resolved Hide resolved

jzuhone added 2 commits August 16, 2019 11:27

Bugfix

c2f90f9

Move this here instead

385a5fa

munkm reviewed Aug 17, 2019

View reviewed changes

jzuhone added 2 commits August 19, 2019 09:23

Make this clearer

4105436

Doc fixes

a116e4a

munkm approved these changes Aug 19, 2019

View reviewed changes

munkm merged commit 04ec81f into yt-project:yt-4.0 Aug 19, 2019

jzuhone deleted the astro_fields branch October 10, 2019 13:44

chummels mentioned this pull request May 19, 2021

H/He phase information not accurate for simple datasets trident-project/trident#165

Closed

jzuhone mentioned this pull request May 20, 2021

Refactoring "default" species fields #3288

Merged

2 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Some updates and additions to astrophysics-related fields #2245

Some updates and additions to astrophysics-related fields #2245

jzuhone commented Apr 17, 2019 •

edited

jzuhone commented Jul 5, 2019

brittonsmith left a comment

jzuhone commented Jul 8, 2019 •

edited

jzuhone commented Jul 8, 2019

brittonsmith commented Jul 9, 2019

AshKelly commented Jul 9, 2019

jzuhone commented Jul 9, 2019

matthewturk commented Jul 12, 2019

jzuhone commented Jul 12, 2019

brittonsmith Jul 16, 2019

jzuhone Jul 18, 2019

brittonsmith Jul 16, 2019

jzuhone Jul 18, 2019

brittonsmith Jul 16, 2019

jzuhone Jul 18, 2019

brittonsmith Jul 16, 2019

brittonsmith Jul 16, 2019

brittonsmith left a comment

brittonsmith Jul 23, 2019

jzuhone commented Aug 8, 2019

jzuhone commented Aug 14, 2019

matthewturk left a comment

matthewturk commented Aug 16, 2019

jzuhone commented Aug 16, 2019

munkm left a comment

munkm Aug 17, 2019

jzuhone Aug 19, 2019

jzuhone commented Aug 19, 2019

munkm left a comment

Some updates and additions to astrophysics-related fields #2245

Some updates and additions to astrophysics-related fields #2245

Conversation

jzuhone commented Apr 17, 2019 • edited

PR Summary

PR Checklist

jzuhone commented Jul 5, 2019

brittonsmith left a comment

Choose a reason for hiding this comment

jzuhone commented Jul 8, 2019 • edited

jzuhone commented Jul 8, 2019

brittonsmith commented Jul 9, 2019

AshKelly commented Jul 9, 2019

jzuhone commented Jul 9, 2019

matthewturk commented Jul 12, 2019

jzuhone commented Jul 12, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

brittonsmith left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jzuhone commented Aug 8, 2019

jzuhone commented Aug 14, 2019

matthewturk left a comment

Choose a reason for hiding this comment

matthewturk commented Aug 16, 2019

jzuhone commented Aug 16, 2019

munkm left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jzuhone commented Aug 19, 2019

munkm left a comment

Choose a reason for hiding this comment

jzuhone commented Apr 17, 2019 •

edited

jzuhone commented Jul 8, 2019 •

edited