[yt-users] potential unit error in cosmology.py

Nathan Goldbaum nathan12343 at gmail.com
Tue Aug 9 05:02:44 PDT 2016


Can one of you file an issue that summarizes the problem so we don't lose
track?

On Tuesday, August 9, 2016, Britton Smith <brittonsmith at gmail.com> wrote:

> Hi Pengfei,
>
> Firstly, Nathan, JohnZ, Matt, Cameron, if you could have a look at the
> last paragraph of this response, I would appreciate it.
>
> Ah, I see now.  I actually get a different answer for the redshift
> interval when I do it by hand.  I get 9.994 instead of 9.979.  Still, our
> answers, differ from the value calculated by the code for the reason you
> describe, essentially that the comoving to proper conversion is different
> between a cosmology object created from a dataset (where a specific
> redshift is defined) and one created by hand.  In the one created from a
> dataset, the comoving and proper frame differ by the factor of (1+z),
> whereas in the one created by hand, they are the same for the reasons I
> described in my previous email.
>
> This is definitely a problem.  We cannot be doing cosmology distance
> calculations in two different unit systems, one where comoving and proper
> are the same and one where they are different.  As strange as this may
> sound, the one that is least correct is using the system where the comoving
> and proper frames are different.  That's because something like a comoving
> radial distance can only be expressed in the comoving frame and has no
> direct translation to the proper frame via some factor of (1 + z).
>
> The solution to this, as I see it, is the following:
> The cosmology object that gets associated with a dataset must have its
> unit_registry overridden to set the comoving and proper frames equal to
> each other.  A better solution would be to even remove the proper frame
> entirely from the cosmology object's unit_registry, though I don't even
> know if that's possible.  I would appreciate the input of someone with more
> intimate knowledge of the unit system, like Nathan, Matt, or John Z.  As a
> quick test, I experimented with the case where we create a cosmology object
> from a dataset by passing a copy of the unit_registry and then modifying
> the comoving frame.  I ran into some difficulties because it looks like,
> for example, Mpc is define independently of pc, and so I would have to
> modify both pccm and Mpcm, and perhaps a number of others that I am unaware
> of.  If anyone has ideas on how best to do this, please say so.
>
> In conclusion, I do think there is an issue that needs to be fixed.  Input
> here would be appreciated.  Thanks, Pengfei, for identifying this issue!
>
> Britton
>
> On Mon, Aug 8, 2016 at 10:25 PM, Pengfei Chen <madcpf at gmail.com
> <javascript:_e(%7B%7D,'cvml','madcpf at gmail.com');>> wrote:
>
>> Hi Britton,
>>
>> Thanks very much for your test! The error happens when I make light rays
>> from one single dataset. I tried to reproduce the error with the dataset
>> under src/yt-hg/tests. The code I used is here:
>>
>> http://paste.yt-project.org/show/6751/
>> After running the code, I did
>>
>> h5ls -d lightray.h5/grid/redshift
>> and found the last redshift was 9.98999107782409. However, the correct
>> value should be ~ 9.979.
>>
>> When I call LightRay this way, it will call the _deltaz_forward function
>> in cosmology_splice.py, which will call the comoving_radial_distance in
>> cosmology.py:
>>
>>         distance2 = self.cosmology.comoving_radial_distance(z2, z)
>> , where distance2 is not in comoving units. Then it will calculate z2:
>>
>>             z2 = ((target_distance - distance2) / m.in_units("Mpccm /
>> h")) + z2
>>
>> . I think the subtraction here is the problem, since it tries to convert
>> distance2 to comoving units, even though the numerical value of distance2
>> is already comoving. So a wrong factor of (1+z) is multiplied.
>>
>> Any comment is appreciated.
>>
>> Thanks again for your help!
>> Pengfei
>>
>> On Mon, Aug 8, 2016 at 7:10 AM, Britton Smith <brittonsmith at gmail.com
>> <javascript:_e(%7B%7D,'cvml','brittonsmith at gmail.com');>> wrote:
>>
>>> Hi Pengfei,
>>>
>>> I'm not able to reproduce what you're seeing, but before I go into that,
>>> I can explain a few things about how the cosmology calculator works.  In
>>> yt, we are only able to define a comoving unit system in relation to a
>>> proper unit system at a specific redshift.  We are unable to define a
>>> comoving unit system that is the base unit system.  This is a problem for
>>> the cosmology calculator, which only calculates comoving quantities.
>>> Therefore, in this case, we are forced to set the comoving and proper unit
>>> systems to be the same.  You can see this for yourself by running the
>>> following script:
>>>
>>> from yt.utilities.cosmology import Cosmology
>>> co = Cosmology()
>>> print (co.comoving_radial_distance(1.92, 2).to("Mpc/h"))
>>> print (co.comoving_radial_distance(1.92, 2).to("Mpccm/h"))
>>>
>>> Both will be the same answer.
>>>
>>> Now, back to the issue you're seeing.  I ran the following script:
>>> http://paste.yt-project.org/show/6741/
>>> to calculate the redshifts dumps needed for a 80 Mpc/h box to go from z
>>> = 2 to z = 1.92 and got the following result:
>>> CosmologyOutputRedshift[0] = 2.000
>>> CosmologyOutputRedshift[1] = 1.926
>>>
>>> This is pretty close to what you expected.  Additionally, if I run the
>>> following script,
>>> http://paste.yt-project.org/show/6743/
>>> with the above redshifts added to the parameter file, I get that dz_max
>>> (the change in redshift when moving 80 Mpccm/h at z = 0)
>>> of 0.0746521607364, which also seems to agree with what you are finding.
>>>
>>> Can you share with me the simulation parameter file and script that
>>> you're using to make your light rays?  I can take a look at that and see if
>>> I still get what I expect.
>>>
>>> Britton
>>>
>>> On Mon, Aug 8, 2016 at 6:13 AM, Pengfei Chen <madcpf at gmail.com
>>> <javascript:_e(%7B%7D,'cvml','madcpf at gmail.com');>> wrote:
>>>
>>>> Hi all,
>>>>
>>>> I was trying to generate light rays using yt 3.4-dev, but found that
>>>> the light rays have wrong redshift intervals. For example, my simulation
>>>> box is 80Mpccm/h per side, and I expect a light ray generated in a data
>>>> dump at z=2 with length =1 simulation unit (i.e. 80Mpccm/h) would have a
>>>> redshift interval from z=2 to z=1.92. However, the light ray generated by
>>>> yt gives spans from z=2 to z=1.97.
>>>>
>>>> As I look into this, I found that the function comoving_radial_distance
>>>> in cosmology.py might return the wrong unit. I think the return value
>>>> should be in comoving units, instead of physical units. To see this
>>>> directly, I made the following change in cosmology_splice.py:
>>>>
>>>> (root) ~/yt-conda/src/yt-hg/yt/analysis_modules/cosmological_observation
>>>> $diff cosmology_splice.py cosmology_splice.py0
>>>> 373d372
>>>> <             print target_distance, distance2, z2
>>>>
>>>> And it shows:
>>>>
>>>> 80.0 Mpccm/h 4.65320708035e+26 cm 1.97387969592 dimensionless
>>>> 80.0 Mpccm/h 1.19455177073e+26 cm 1.97343368065 dimensionless
>>>> 80.0 Mpccm/h 1.21507519065e+26 cm 1.97342592996 dimensionless
>>>>
>>>> Then I made the following change in cosmology.py:
>>>>
>>>> (root) ~/yt-conda/src/yt-hg/yt/utilities $diff cosmology.py
>>>> cosmology.py0
>>>> 111,112c111,112
>>>> <         return self.quan((self.hubble_distance() *
>>>> <                 trapzint(self.inverse_expansion_factor, z_i,
>>>> z_f)).value, 'cmcm')
>>>> ---
>>>> >         return (self.hubble_distance() *
>>>> >                 trapzint(self.inverse_expansion_factor, z_i,
>>>> z_f)).in_base(self.unit_system)
>>>>
>>>>
>>>> Then I get:
>>>> 80.0 Mpccm/h 4.65320708035e+26 cmcm 1.92163908776 dimensionless
>>>> 80.0 Mpccm/h 3.62771515661e+26 cmcm 1.92124702027 dimensionless
>>>>
>>>> With the change of unit from cm to cmcm, the light rays have the right
>>>> span.
>>>>
>>>> Even though this solves my problem, I am not sure if similar problem
>>>> still exists. For example, instead of making change in
>>>> comoving_radial_distance, we might need to change hubble_distance into
>>>> comoving units. Hopefully someone familiar with yt unit system could check
>>>> this.
>>>>
>>>> Thanks,
>>>> Pengfei
>>>>
>>>> _______________________________________________
>>>> yt-users mailing list
>>>> yt-users at lists.spacepope.org
>>>> <javascript:_e(%7B%7D,'cvml','yt-users at lists.spacepope.org');>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>
>>>>
>>>
>>> _______________________________________________
>>> yt-users mailing list
>>> yt-users at lists.spacepope.org
>>> <javascript:_e(%7B%7D,'cvml','yt-users at lists.spacepope.org');>
>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>
>>>
>>
>> _______________________________________________
>> yt-users mailing list
>> yt-users at lists.spacepope.org
>> <javascript:_e(%7B%7D,'cvml','yt-users at lists.spacepope.org');>
>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20160809/fd13ce1b/attachment.htm>


More information about the yt-users mailing list