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

Pengfei Chen madcpf at gmail.com
Tue Aug 9 08:12:06 PDT 2016


Hi Nathan,

Thanks for your reminding. I just opened an issue:
https://bitbucket.org/yt_analysis/yt/issues/1258/potential-unit-error-in-cosmologypy

Thanks,
Pengfei

On Tue, Aug 9, 2016 at 5:02 AM, Nathan Goldbaum <nathan12343 at gmail.com>
wrote:

> 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> 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>
>>> 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> 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
>>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>>
>>>>>
>>>>
>>>> _______________________________________________
>>>> yt-users mailing list
>>>> 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
>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>
>>>
>>
> _______________________________________________
> yt-users mailing list
> 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/ff0b91c7/attachment-0001.html>


More information about the yt-users mailing list