[yt-dev] [continued] RAMSES units bug in yt-3.0

Matthew Turk matthewturk at gmail.com
Tue May 13 05:52:06 PDT 2014 

On Tue, May 13, 2014 at 8:48 AM, Sam Geen <samgeen at astro.ox.ac.uk> wrote:
> Ah, very good, thanks for the link to that old PR! I was worried I was going
> crazy... in any case, whichever is the easiest way to make the
> position/length units consistent.

I change the particle positions in the ingestion phase, so at no point
should they appear incorrect.  I'd really like to figure out the
boxlen != 1.0 ...

>
> By the way, I took a quick look at reading the grav_ files recently,
> although the output seems kind of strange, so I haven't issued a PR yet
> until I understand what it's doing (plus I still don't 100% understand the
> units). I suspect it might be reading in the wrong order or something. It
> also complicates the code, so it might be worth refactoring a bit to make a
> universal "AMR-like" data_structures/io interface for the files called
> hydro_, grav_, rt_, etc, rather than what I did which was to scatter a bunch
> of conditionals into the existing functions.

That could be -- I think there was some discussion that they might
also include ghost zones?

I've also run some checks on cosmology runs.  Here's what I've found.
With no changes to the current tip of yt-3.0, I have run this script:

--
ds = yt.load("output_00010/info_00010.txt")
dd = ds.all_data()

# Let's find the total mass
pmass = dd["particle_mass"].sum().in_units("Msun")
cmass = dd["cell_mass"].sum().in_units("Msun")
rho = ds.cosmology.critical_density(ds.current_redshift)
crit_mass = (rho * dd["cell_volume"]).sum().in_units("Msun")
print pmass
print cmass
print pmass / (pmass + cmass)
print ((pmass + cmass) / crit_mass)
--

The results are good.  I get roughly 0.84 for pmass / (pmass + cmass),
and I get 0.998 for the final check, which seems within the parameters
of the simulation.  So the cosmological units should be correct, and
more to the point, they match up with what I'm expecting having read
over the code carefully.

So it looks to me like there is in fact a boxlen missing, but that as
it stands, the particle and density units are correct for boxlen =
1.0.  Nick's message is worrisome because the factor isn't exactly 1e3
between his max density.  I have a boxlen != 1.0 dataset, but I have
no reference values to compare it again.  Any chance there's one with
*reference* results for total mass in particles, total mass in gas,
etc?

-Matt

>
>
> On 13/05/2014 14:21, Matthew Turk wrote:
>>
>> Hi Sam,
>>
>> So, here's an old pull request I found:
>>
>>
>> http://hg.yt-project.org/yt-3.0/pull-request/62/adding-boxlen-to-ramses-units-for-mass-and
>>
>> and then there's the units.f90 from RAMSES itself:
>>
>>
>> https://bitbucket.org/rteyssie/ramses/src/81ab29b8a405e7ccf6bf30fd67582b4edacddd6e/trunk/ramses/amr/units.f90?at=master
>>
>> I'm going to take care of this, update the PR, and ping you both.
>> (And maybe we should move to yt-dev? :)
>>
>> -Matt
>>
>> On Tue, May 13, 2014 at 8:17 AM, Sam Geen <samgeen at astro.ox.ac.uk> wrote:
>>>
>>> Huh, that difference in the max density is very strange. I can't think
>>> why
>>> you'd have that difference; it's not an obvious multiple.
>>>
>>> Actually, I'm looking at the code and it's possible that the particle
>>> positions are scaled from 0 to boxlen and the grid positions from 0 to 1,
>>> though I'd need to open up the raw values in a RAMSES file to confirm
>>> this.
>>> If this is true, I suppose one "fix" could be to hard-multiply the cell
>>> sizes by pf["boxlen"] when you load them so that the length units are
>>> consistent for all data. I *think* that the older version of yt-3.0 gave
>>> correct answers, though, so that should at least be a way to calibrate
>>> the
>>> new version.
>>>
>>>
>>> On 13/05/2014 14:01, nick moeckel wrote:
>>>
>>> I had a chance to test out an older and a newer version of yt-3 on a
>>> small
>>> dataset that has boxlen=10. This seems to confirm a boxlen**3 factor
>>> somewhere, although there's a further difference in the max density.
>>>
>>> more recent version (hg id -i gives f4838a2165c0):
>>>
>>> after dd = ds.all_data():
>>>
>>> In [11]: dd.quantities.extrema('density')
>>>
>>> Out[11]: (2.00739832795e-28 g/cm**3, 4.54771338163e-24 g/cm**3)
>>>
>>>
>>> older version (hg id -i gives 3e8b733c9ee9):
>>>
>>> after dd = ds.h.all_data():
>>>
>>> In [10]: dd.quantities.extrema('Density')
>>>
>>> Out[10]: [(2.0073983279492632e-25, 4.0079147772594363e-21)]
>>>
>>>
>>>
>>> On Fri, May 9, 2014 at 5:50 PM, Sam Geen <samgeen at astro.ox.ac.uk> wrote:
>>>>
>>>> Yep, that looks like it should work. I'll try to run it on some particle
>>>> data when I get the time, but like I said I'm 99% sure the mass units
>>>> should
>>>> be identical for both grid hydro and particle data.
>>>>
>>>>
>>>> On 09/05/14 17:46, Matthew Turk wrote:
>>>>>
>>>>> Hi Sam,
>>>>>
>>>>> Okay, I've looked over a bit, and I think the correct change would be:
>>>>>
>>>>> mass_unit = rho_u * length_unit**3.
>>>>>
>>>>> That should include the correct length unit, and I think will reduce
>>>>> "density" back to the "unit_d" that's in the parameter file.  If this
>>>>> looks okay to you, I will push it, but I really do want to make sure
>>>>> the particle masses are correct.  Can Nick or Romain provide a bit of
>>>>> guidance here?
>>>>>
>>>>> -Matt
>>>>>
>>>>> On Fri, May 9, 2014 at 11:25 AM, Sam Geen <samgeen at astro.ox.ac.uk>
>>>>> wrote:
>>>>>>
>>>>>> OK, interesting. In theory RAMSES should have identical units for both
>>>>>> particles and gas. I can hunt down a run with particles to test if you
>>>>>> like.
>>>>>>
>>>>>> Thanks!
>>>>>> Sam
>>>>>>
>>>>>>
>>>>>> On 09/05/14 17:22, Matthew Turk wrote:
>>>>>>>
>>>>>>> Hi Sam,
>>>>>>>
>>>>>>> On reflection, I think this might be related to getting the
>>>>>>> *particle*
>>>>>>> masses correct.  I will take a look at it as soon as I can.
>>>>>>>
>>>>>>> -Matt
>>>>>>>
>>>>>>> On Fri, May 9, 2014 at 11:00 AM, Sam Geen <samgeen at astro.ox.ac.uk>
>>>>>>> wrote:
>>>>>>>>
>>>>>>>> Hmm, that looks like it should be "mass_unit = rho_u *
>>>>>>>> length_unit**3"
>>>>>>>> in
>>>>>>>> line 492. You're right that it mentions the boxlength issue, though.
>>>>>>>>
>>>>>>>>
>>>>>>>> On 09/05/14 16:50, Matthew Turk wrote:
>>>>>>>>>
>>>>>>>>> Hi Sam,
>>>>>>>>>
>>>>>>>>> Okay, sounds good.  Looking at how code unit attributes are set up:
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> https://bitbucket.org/yt_analysis/yt/src/a14a150c7c81850df81346162bdaff271e77eb50/yt/frontends/ramses/data_structures.py?at=yt-3.0#cl-482
>>>>>>>>>
>>>>>>>>> suggests to me that length_unit takes into account boxlen, and
>>>>>>>>> mass_unit does not.  The comments have some indication why this
>>>>>>>>> might
>>>>>>>>> be.
>>>>>>>>>
>>>>>>>>> -Matt
>>>>>>>>>
>>>>>>>>> On Fri, May 9, 2014 at 10:46 AM, Sam Geen <samgeen at astro.ox.ac.uk>
>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>> Yep; both this and dd["density"] give cgs values that are too
>>>>>>>>>> small
>>>>>>>>>> by
>>>>>>>>>> roughly a factor of boxlen**3.
>>>>>>>>>>
>>>>>>>>>> One other thing I need to try is to make sure I'm using the very
>>>>>>>>>> latest
>>>>>>>>>> version of YT; I've been playing around with the Ramses frontend
>>>>>>>>>> so
>>>>>>>>>> it's
>>>>>>>>>> possible my version is somehow out of sync. Will let you know if
>>>>>>>>>> this
>>>>>>>>>> fixes
>>>>>>>>>> things.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> On 09/05/14 16:42, Matthew Turk wrote:
>>>>>>>>>>>
>>>>>>>>>>> Hi Sam,
>>>>>>>>>>>
>>>>>>>>>>> Can you verify the units are in fact incorrect in *cgs*?
>>>>>>>>>>> Something
>>>>>>>>>>> like this would work:
>>>>>>>>>>>
>>>>>>>>>>> ds = load(...)
>>>>>>>>>>> dd = ds.all_data()
>>>>>>>>>>> print dd.quantities.total_quantitiy("cell_mass").in_cgs()
>>>>>>>>>>>
>>>>>>>>>>> Your second message makes me wonder if there's just a slipup in
>>>>>>>>>>> how
>>>>>>>>>>> the units are returned.
>>>>>>>>>>>
>>>>>>>>>>> -Matt
>>>>>>>>>>>
>>>>>>>>>>> On Fri, May 9, 2014 at 10:37 AM, Sam Geen
>>>>>>>>>>> <samgeen at astro.ox.ac.uk>
>>>>>>>>>>> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>> Sorry for the spam; a second bug I've seen is that the density
>>>>>>>>>>>> and
>>>>>>>>>>>> pressure
>>>>>>>>>>>> unit labels on figures appears to be broken; it seems to print a
>>>>>>>>>>>> latex-mangled code names for the units rather than the cgs name
>>>>>>>>>>>> of
>>>>>>>>>>>> the
>>>>>>>>>>>> units; see attached example. Temperature is fine.
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> On 09/05/14 16:28, Sam Geen wrote:
>>>>>>>>>>>>>
>>>>>>>>>>>>> Hi,
>>>>>>>>>>>>>
>>>>>>>>>>>>> Just to say that in the latest version of yt-3.0 (i.e. since
>>>>>>>>>>>>> various
>>>>>>>>>>>>> fields were renamed or re-implemented), I've found a bug in the
>>>>>>>>>>>>> implementation of "cell_mass", which is giving results that are
>>>>>>>>>>>>> too
>>>>>>>>>>>>> low
>>>>>>>>>>>>> in
>>>>>>>>>>>>> my runs; I believe the issue is that it's missing a factor of
>>>>>>>>>>>>> pf["boxlen"]**3 (which is of course only a problem if boxlen is
>>>>>>>>>>>>> not
>>>>>>>>>>>>> 1).
>>>>>>>>>>>>> The
>>>>>>>>>>>>> previous "CellMass_Msun" worked fine. If I get time I might
>>>>>>>>>>>>> take
>>>>>>>>>>>>> a
>>>>>>>>>>>>> look
>>>>>>>>>>>>> and
>>>>>>>>>>>>> issue a pull request, but otherwise I'm just flagging this in
>>>>>>>>>>>>> case
>>>>>>>>>>>>> someone
>>>>>>>>>>>>> else runs into problems; you can just manually multiply the
>>>>>>>>>>>>> result
>>>>>>>>>>>>> by
>>>>>>>>>>>>> pf["boxlen"]**3 until it's fixed.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Cheers,
>>>>>>>>>>>>>
>>>>>>>>>>>>> Sam
>>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>>> 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
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> 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
>>>>>
>>>>> _______________________________________________
>>>>> 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
>>>
>> _______________________________________________
>> yt-users mailing list
>> yt-users at lists.spacepope.org
>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
>
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org

More information about the yt-dev mailing list