[yt-users] Problem with Units in a Derived Field

Jared Coughlin Jared.W.Coughlin.29 at nd.edu
Wed Aug 26 08:34:17 PDT 2015


My ultimate goal is to make an absorption spectrum, so I was trying to add
an HI number density field.  I don't have any experience working with grid
based codes, only sph, so I was assuming that I was using the particle
positions.
-Jared

On Wed, Aug 26, 2015 at 11:25 AM, Nathan Goldbaum <nathan12343 at gmail.com>
wrote:

> So the issue is that you're combining particle fields, which are defined
> at the location of your SPH particles, with fields that are defined on the
> octree mesh.
>
> data['Gas', 'InternalEnergy']
>
> ^ is a particle field defined at the locations of your SPH particles
> (particles of type 'Gas').
>
> data['gas', 'density']
>
> ^ is a mesh field defined at the cell centers of the octree mesh yt
> constructs in memory. It's specifically generated an SPH smoothing
> operation.
>
> The reasons the shapes of the arrays you're getting back are different is
> that the number of leaf octree cells is not the same as the number of Gas
> particles.
>
> I'm not sure if what you're trying to do is supposed to use the mesh
> locations or the particle locations. The correct way to modify your script
> would be different in the two cases.
>
> On Wed, Aug 26, 2015 at 10:17 AM, Jared Coughlin <
> Jared.W.Coughlin.29 at nd.edu> wrote:
>
>> Sorry about that! Here's the script:
>>
>> http://paste.yt-project.org/show/5842/
>>
>> And the traceback:
>>
>> http://paste.yt-project.org/show/5843/
>>
>> Thanks!
>> -Jared
>>
>> On Tue, Aug 25, 2015 at 7:41 PM, Nathan Goldbaum <nathan12343 at gmail.com>
>> wrote:
>>
>>>
>>>
>>> On Tuesday, August 25, 2015, Jared Coughlin <Jared.W.Coughlin.29 at nd.edu>
>>> wrote:
>>>
>>>> Hello! I'm really sorry for the barrage of questions; hopefully this
>>>> will be the last one!
>>>> I now get the error:  operands could not be broadcast together with
>>>> shapes (2072685,) (2153824,)
>>>> on line: fHI = (recomb * ne) / (gamma_HI + (gamma_c * ne))
>>>>
>>>> I printed the shapes of all of those arrays, and this is what I get:
>>>> nH.shape
>>>> (16, 16, 16)
>>>> ne.shape
>>>> (16, 16, 16)
>>>> gamma_c.shape
>>>> (1,)
>>>> recomb.shape
>>>> (1,)
>>>> nH.shape
>>>> (2153824,)
>>>> ne.shape
>>>> (2153824,)
>>>> gamma_c.shape
>>>> (2072685,)
>>>> recomb.shape
>>>> (2072685,)
>>>>
>>>> Firstly, I'm a little confused as to why their dimensions change, but
>>>> what really matters are the bottom ones.  I know how many gas particles are
>>>> in my gadget simulation: 2072685.  All of the arrays were created from
>>>> either data["Gas", "InternalEnergy"] or data["gas", "density"], which
>>>> should, I would imagine, have the same dimensions: the number of gas
>>>> particles.  I was just wondering if anyone had seen something like this
>>>> before? The function in question was posted to a pastebin in an earlier
>>>> email on this thread, and I have since made all of the arrays YTArrays, as
>>>> suggested above. Thanks!
>>>> -Jared
>>>>
>>>>
>>> Hey Jared,
>>>
>>> I looked above but can't seem to find the correct link for the script
>>> you're using. Can you pastebin your updated script and the full error
>>> tracrback you're seeing now?
>>>
>>> Nathan
>>>
>>>
>>>
>>>>
>>>> On Tue, Aug 25, 2015 at 5:58 PM, Jared Coughlin <
>>>> Jared.W.Coughlin.29 at nd.edu> wrote:
>>>>
>>>>> Thanks, I'll try that!
>>>>> -Jared
>>>>>
>>>>> On Tue, Aug 25, 2015 at 5:55 PM, Andrew Myers <atmyers2 at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Hi Jared,
>>>>>>
>>>>>> You probably don't want to loop over the entire array in pure python;
>>>>>> that could be very slow for large datasets. One thing you can do is use
>>>>>> slicing to select the elements of the array that meet your criterion:
>>>>>>
>>>>>> indices = data['gas','InternalEnergy'] <= 1e4
>>>>>> data['gas', 'InternalEnergy'][indices] = *something*
>>>>>>
>>>>>> that would set all the elements where the InternalEnergy field is
>>>>>> less than or equal to 1e4 to *something*. Usually you can express most
>>>>>> operations you'd want to do this way without explicitly looping over the
>>>>>> elements.
>>>>>>
>>>>>> -Andrew
>>>>>>
>>>>>> On Tue, Aug 25, 2015 at 2:44 PM, Jared Coughlin <
>>>>>> Jared.W.Coughlin.29 at nd.edu> wrote:
>>>>>>
>>>>>>> That seemed to have worked, thanks!
>>>>>>>
>>>>>>> I just have another quick question, if you wouldn't mind? In my
>>>>>>> function to generate my new field, there's an if statement that looks at
>>>>>>> the value of the temperature: if data['gas','InternalEnergy'] <= 1e4:
>>>>>>> However, I get the error: "The truth value of an array with more
>>>>>>> than one element is ambiguous. Use a.any() or a.all()." I wasn't sure how
>>>>>>> to loop over the elements of the array to access a specific element, or if
>>>>>>> that's even possible? Using a.any() or a.all() don't seem to fit the bill
>>>>>>> for this particular problem, either. Thanks!
>>>>>>> -Jared
>>>>>>>
>>>>>>> On Tue, Aug 25, 2015 at 3:19 PM, Nathan Goldbaum <
>>>>>>> nathan12343 at gmail.com> wrote:
>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> On Tue, Aug 25, 2015 at 2:17 PM, Jared Coughlin <
>>>>>>>> Jared.W.Coughlin.29 at nd.edu> wrote:
>>>>>>>>
>>>>>>>>> Great, thanks! I tried that, and it gave me the error: "YTQuantity
>>>>>>>>> instances must be scalars." I'm not sure what to do about that, since all
>>>>>>>>> of these quantities depend on either temperature or density, both of which
>>>>>>>>> are arrays. I tried looking at the doc page for YTQuantity, but I didn't
>>>>>>>>> see anything that looked like it could solve this problem, though it's
>>>>>>>>> possible I missed it.
>>>>>>>>>
>>>>>>>>
>>>>>>>> You need to use YTArray instead. YTQuantity is a subclass of
>>>>>>>> YTArray that only holds scalars. YTArray can hold arrays of data with units.
>>>>>>>>
>>>>>>>>
>>>>>>>>> -Jared
>>>>>>>>>
>>>>>>>>> On Mon, Aug 24, 2015 at 8:13 PM, Nathan Goldbaum <
>>>>>>>>> nathan12343 at gmail.com> wrote:
>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> On Monday, August 24, 2015, Cameron Hummels <chummels at gmail.com>
>>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>>> Hi Jared,
>>>>>>>>>>>
>>>>>>>>>>> The problem is one of units.  In order to make common arithmetic
>>>>>>>>>>> work like addition, yt requires both arguments going into the addition
>>>>>>>>>>> operator to have the same units.
>>>>>>>>>>>
>>>>>>>>>>> In your line,
>>>>>>>>>>>
>>>>>>>>>>> fHI = (recomb * ne) / (gamma_HI + (gamma_c *ne))
>>>>>>>>>>>
>>>>>>>>>>> you're adding two things (gamma_HI) and (gamma_c*ne) that have
>>>>>>>>>>> different units--in this case, something with (code_mass/code_length**3)
>>>>>>>>>>> units and something with no units defined as (1).
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Just a quick clarification: they need not have the same units,
>>>>>>>>>> just the same dimensions.  You can add two quantities with units of grams
>>>>>>>>>> and solar masses (for example) but not grams and liters.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> In order to fix this, you can re-assign the appropriate units to
>>>>>>>>>>> different arrays or quantities with the YTArray and YTQuantity classes.  In
>>>>>>>>>>> this case, make sure that the two arguments going into the addition have
>>>>>>>>>>> the same units.  If you want to recast "ne" to have number density units,
>>>>>>>>>>> you can do this:  from yt.units.yt_array import YTQuantity; ne =
>>>>>>>>>>> YTQuantity(ne, 'cm**-3').  For more info on units, check out
>>>>>>>>>>> http://yt-project.org/docs/dev/analyzing/units/index.html
>>>>>>>>>>>
>>>>>>>>>>> I hope this helps!
>>>>>>>>>>>
>>>>>>>>>>> Cameron
>>>>>>>>>>>
>>>>>>>>>>> On Mon, Aug 24, 2015 at 4:12 PM, Jared Coughlin <
>>>>>>>>>>> Jared.W.Coughlin.29 at nd.edu> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> Hello! I have a gadget snapshot for which the standard internal
>>>>>>>>>>>> energy block has been replaced by one containing the temperature.  I'm
>>>>>>>>>>>> trying to calculate a derived field for the neutral hydrogen number
>>>>>>>>>>>> density,  but I'm getting an error when I try to do ds.add_field():
>>>>>>>>>>>>
>>>>>>>>>>>> http://paste.yt-project.org/show/5833/
>>>>>>>>>>>>
>>>>>>>>>>>> However, I get an error about being unable to add quantities
>>>>>>>>>>>> with differing units:
>>>>>>>>>>>>
>>>>>>>>>>>> http://paste.yt-project.org/show/5834/
>>>>>>>>>>>>
>>>>>>>>>>>> The docs say not to do any unit conversion because that is
>>>>>>>>>>>> apparently taken care of behind the scenes, so I didn't.  I was just
>>>>>>>>>>>> wondering if there was a way around this? Thanks!
>>>>>>>>>>>> -Jared
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>> yt-users mailing list
>>>>>>>>>>>> yt-users at lists.spacepope.org
>>>>>>>>>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> --
>>>>>>>>>>> Cameron Hummels
>>>>>>>>>>> NSF Postdoctoral Fellow
>>>>>>>>>>> Department of Astronomy
>>>>>>>>>>> California Institute of Technology
>>>>>>>>>>> http://chummels.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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20150826/114436c6/attachment.html>


More information about the yt-users mailing list