Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[enhancement] Add relativistic component to cosmology calculations #1714

Merged
merged 28 commits into from Feb 11, 2019

Conversation

brittonsmith
Copy link
Member

@brittonsmith brittonsmith commented Mar 7, 2018

Update: I've finally come back to this and it should now be in a finished state.

This adds the relativistic component of the energy density of the universe, omega_radiation, to the cosmology calculations. Since there is no analytic solution to the Friedmann equations when omega_radiation > 0, we now numerically integrate using the trapezoid rule and then linearly interpolate in log(a) vs. log(t).

To get z_from_t, we reverse that table and interpolate. In this case, I also iteratively adjust the bounds of the log(a) table to be outside of the calculated log(a) values. This makes sure the integration is accurate and the ends of the table. This is most important at the high redshift end.

It's also important to note that the previous z_from_t and t_from_z calculations did not support all cosmologies, whereas this does.

In response to @jzuhone's comment about the hubble_time calculation, I've switched it to the correct formulation. @jzuhone was definitely right all along about that.

Differences of redshift/time conversions (t_from_z and z_from_t) are maximally about 1e-5. It's not nothing, but I think quite acceptable.

This is now supported in the Enzo code, so that frontend is updated to include it. I have added it as a standard dataset attribute for cosmology simulations.

@jzuhone
Copy link
Contributor

jzuhone commented Mar 7, 2018

So this is a bit tangential, but I hadn't really looked at the hubble_time method until now, and now I realize that what it actually computes is the age of the universe. When I hear "hubble time" I think of 1/H_0.

Would it make sense to address this naming issue in this PR as well? Maybe rename the hubble_time method to something about the age of the universe and then let hubble_time point to it with a deprecation warning?

Or am I missing something here entirely?

@brittonsmith
Copy link
Member Author

I will use "hubble time" in a redshift-dependent way in speech, as in "the hubble time at redshift 20", but maybe that's just me. Does anyone else do this?

@jzuhone
Copy link
Contributor

jzuhone commented Mar 7, 2018

I've always seen it as "Hubble time" = t_H = 1/H_0. For a flat matter-only universe, the age of the universe is t = 2/3t_H = 2/(3H_0). But the hubble_time method is clearly calculating the age and not the "Hubble time".

@jzuhone
Copy link
Contributor

jzuhone commented Mar 7, 2018

Obviously "Hubble time" is redshift-dependent since H is a parameter which depends on z, but the method itself seems to compute not this but the age.

@brittonsmith
Copy link
Member Author

Alright, what do you suggest then? Should hubble_time be deprecated in favor of something like universe_age?

@jzuhone
Copy link
Contributor

jzuhone commented Mar 8, 2018

I think that works, but open to suggestions (also open to challenges from anyone who thinks I am completely full of it)

@ngoldbaum
Copy link
Member

ngoldbaum commented Mar 8, 2018

Is it useful to also have a helper that returns 1/H? Not sure what that would be called if hubble_time is taken.

@brittonsmith
Copy link
Member Author

I'm willing to say @jzuhone is correct, in which case 1/H should be hubble_time, but that does break things. I suppose we already have hubble_parameter, which returns H(z), so maybe it's not necessary and we should just deprecate hubble_time.

@brittonsmith brittonsmith changed the title [enhancement] Add relativistic component to cosmology calculations [wip, enhancement] Add relativistic component to cosmology calculations Apr 10, 2018
@brittonsmith brittonsmith changed the title [wip, enhancement] Add relativistic component to cosmology calculations [enhancement] Add relativistic component to cosmology calculations Feb 4, 2019
@brittonsmith
Copy link
Member Author

I've finally gotten around to finishing this off. It should be ready for review now.

Copy link
Member

@ngoldbaum ngoldbaum left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't reviewed the math in the new a_from_t and t_from_a functions.

Do you think the existing tests for the cosmology calculator are sufficient to be confident about the changes here? It looks like only hubble_time,t_from_z, and z_from_t are explicitly tested. It also looks like hubble_time and t_from_z are tested to be identical in the default cosmology, which I guess is still the case.

Still would be nice to get some more tests, particularly for the new functionality you've added.

Do you have any Enzo datasets with nonzero omega_radiation we could add to yt-project.org/data (if we don't have those already). That would make it easier to debug any issues people report in the future about this. Could also use such a dataset for a test of this functionality.

@@ -319,9 +323,9 @@ def lookback_time(self, z_i, z_f):
return (trapzint(self.age_integrand, z_i, z_f) / \
self.hubble_constant).in_base(self.unit_system)

def hubble_time(self, z, z_inf=1e6):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I realize we talked about this already, but would you be ok with keeping the old calculation but deprecating the name hubble_time. You could then add a new function which returns one over the hubble constant. I worry about users who will now get different results from this function in the next yt version with no warning. If you feel that's not a concern worth worrying about I'd love to hear why you think so.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, I can do that. I think the counterargument is that this could be seen as constituting a bugfix, but since the correct behavior is basically trivial, I can just include that in the deprecation warning.

@@ -388,10 +384,11 @@ def hubble_parameter(self, z):
>>> print(co.hubble_parameter(1.0).in_units("km/s/Mpc"))

"""
return self.hubble_constant.in_base(self.unit_system) * self.expansion_factor(z)
return self.hubble_constant.in_base(self.unit_system) * \
self.expansion_factor(z)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor nit, can you format this using parens to avoid the continuation character:

return (self.hubble_constant.in_base(self.unit_system) *
        self.expansion_factor(z))

@brittonsmith
Copy link
Member Author

@ngoldbaum, I'm happy to add some tests of each function if I'm allowed to just compare with some hard-coded answers. I can also do some comparisons with the old analytic functions from Enzo for the supported cosmologies. I'll also see what I can do about some Enzo data.

@yt-project yt-project deleted a comment from codecov bot Feb 6, 2019
@brittonsmith
Copy link
Member Author

I've issues a PR to the website repo to add an Enzo dataset with non-zero omega_radiation. I've also added some small tests using this dataset here.

@ngoldbaum
Copy link
Member

@yt-fido test this please

@brittonsmith
Copy link
Member Author

These test failures are real. I'll look into them.

@brittonsmith
Copy link
Member Author

The perils of giving keyword arguments as positional arguments...

@brittonsmith
Copy link
Member Author

The rest of the test failures likely result from minor changes in cosmology calculations. I'll investigate and update when I have the chance.

@ngoldbaum
Copy link
Member

It looks like four of the cookbook recipes are raising errors on python2.

I think bumping the answer numbers is warranted if the amswers are sensitive to that and the changes in answers are small.

@brittonsmith
Copy link
Member Author

@yt-fido test this please

@brittonsmith
Copy link
Member Author

brittonsmith commented Feb 10, 2019

I traced the answer test failures down to a change that I never intended to keep. That's fixed and they now all pass without having to update answers. The only remaining failures are due to pyaml not being installed for the minimal dependencies tests. Would it be ok to add that to those builds?

Edit: I meant pyyaml, not pyaml. They are different.

@ngoldbaum
Copy link
Member

So the whole point of the minimal dependencies tests is that the test suite runs without errors with a minimal set of dependencies installed. Can you make it so the tests that depend on pyyaml only run if it’s importable?

@ngoldbaum
Copy link
Member

But definitely add it to test_requirements.txt.

@ngoldbaum
Copy link
Member

Awesome, thanks so much for sticking with it on this one! I'm liking all the new tests :)

@ngoldbaum ngoldbaum merged commit 7f0b1de into yt-project:master Feb 11, 2019
@brittonsmith brittonsmith deleted the omegarad branch February 11, 2019 20:07
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants