[yt-users] units of cic mesh output from covering_grid

[Nathan Goldbaum]

On Tue, Jul 1, 2014 at 4:31 PM, Nathan Goldbaum <nathan12343 at gmail.com>
wrote:

> Hi Brendan,
>
> Not sure exactly what's happening.  In general it's difficult to remotely
> debug something like this without a simple test case that one of us can run
> locally.  Can you perhaps make a test that illustrates the issue using one
> of the datasets on yt-project.org/data?
>
>
> On Tue, Jul 1, 2014 at 3:05 PM, Brendan Griffen <
> brendan.f.griffen at gmail.com> wrote:
>
>> I have made a mesh using covering_grid:
>>
>> pf = load(basepath+ext+'.0.hdf5', unit_base = unit_base, bounding_box =
>> bbox)
>>
>> unit_base = {'UnitLength_in_cm'         : 3.08568e+21,
>>              'UnitMass_in_g'            :   1.989e+43,
>>              'UnitVelocity_in_cm_per_s' :      100000}
>>
>>
> FWIW, you shouldn't need to set the units or the bbox by hand, yt can read
> this data from the hdf5 file.
>

Hmm, looks like I was a little hasty earlier.  This statement is based on
the OWLS data format, which is a little more descriptive than the default
gadget hdf5 data format.  In particular, units and bounding box information
are written to disk along with the particle data in the hdf5 files.

If your dataset acts like the public Gadget distribution, you will need to
provide both the bbox and the unit information.


>
>
>>  In [3]: bbox
>>
>> Out[3]:
>>
>> array([[   0.,  100.],
>>
>>        [   0.,  100.],
>>
>>        [   0.,  100.]])
>>
>
This might be the issue, depending on what units the bbox is interpreted to
be in. For what it's worth, the OWLS dataset on yt-project.org/data seems
to be a 25 comoving megaparsec / h box:

$ yt load snapshot_033/snap_033.0.hdf5
IPython 2.1.0 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

In [1]: ds = yt.load('snapshot_033/snap_033.0.hdf5')

In [2]: ds.domain_left_edge
Out[1]: YTArray([ 0.,  0.,  0.]) code_length

In [3]: ds.domain_right_edge
Out[2]: YTArray([ 25.,  25.,  25.]) code_length

In [4]: ds.length_unit
Out[3]: 3.085678e+24 cmcm/h



>  ndim = 512
>>
>> level = int(math.log(ndim,2))
>>
>> cg = pf.covering_grid(level=level,
>> left_edge=[0,0,0],dims=[ndim,ndim,ndim])
>>
>> arr = cg['deposit', 'all_cic']
>>
>> When I do arr.mean() etc. I get the output:
>>
>>  --> Doing dimension: 512
>>
>>         == Mass density ==
>>
>>    --- Max:  2.367e-13 [g/cm^-3]
>>
>>    --- Min:  0.000e+00 [g/cm^-3]
>>
>>    --- Mean: 1.764e-21 [g/cm^-3]
>>
>>         == Number density ==
>>
>>    --- Max:  1.158e+11 [cm^-3]
>>
>>    --- Min:  0.000e+00 [cm^-3]
>>
>>    --- Mean: 8.630e+02 [cm^-3]
>>
>> Written out mesh:
>> /bigbang/data/bgriffen/c2ray/cicfiles/parent/512/density/z0.000.dat
>>
>> I thought it returns the mesh in g/cm^3. Compared to just using the
>> particles number, box width and particle mass I get (mass & number density):
>>
>> orig_mass =
>> (head.massarr[1]*10**10*head.nall[1]/head.hubble)/(head.boxwidth/head.hubble)**3
>>
>> # just total_mass/total_volume where "head" => header info of HDF5 gadget
>> snapshot
>>
>> In [29]: orig_mass*MsolMpc3_to_gcm3
>>
>> Out[29]: 2.701e-30
>>
>> In [30]: orig_mass*MsolMpc3_to_gcm3*coeff
>>
>> Out[30]: 1.321e-06
>>
>> These seem reasonable but the values back from yt don't. Can anyone
>> confirm the yt units of the mesh returned? Does it take care of hubble
>> parameter as read from block? Though this would only account for some minor
>> discrepancy - it is orders of magnitude off at this stage.
>>
> yt should detect your dataset as being cosmological and apply the
appropriate scaling factors.

That said I don't have a cosmological gadget hdf5 dataset to test with, so
I can't easily do an apples-to-apples comparison with what you're trying to
accomplish.

My best advice would be to look at the code in the Gadget frontend using a
debugger and try to figure out what's going wrong. We haven't had a large
number of users for the gadget frontend yet so it's likely there's a lot of
room for improvement. Patches are more than welcome.

If anyone happens to have a reasonably sized (< 1 GB) cosomological gadget
hdf5 dataset, it would be an excellent addition to yt-project.org/data.


>  Thanks.
>>
>> Brendan
>>
>> _______________________________________________
>> 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/20140701/8db7c0ba/attachment.htm>