[yt-users] Setting Units

[John Regan]

Kacper/Michael/Chris,

Thanks a million lads - that makes perfect sense! I hadn't realised the
zeros_like function existed - nice one. And thanks Kacper for debugging the
indexing statement that works a charm!

Cheers lads.

John

On Fri, Sep 26, 2014 at 5:53 PM, Chris Malone <chris.m.malone at gmail.com>
wrote:

> Hi John,
>
> Following up on Kacper's comment, which is what I was going to suggest :):
>
> The 'data' passed to a field is actually a chunk of data, and the chunking
> method depends on how it was implemented in the frontend.  That said, the
> derived field should be somewhat agnostic to 'data's size/shape.  In your
> case, you get around this by saying you want k17rate to be
> np.zeros_like(T), which means to have the same shape as T.
>
> On Fri, Sep 26, 2014 at 8:38 AM, John Regan <johnanthonyregan at gmail.com>
> wrote:
>
>> HI Britton,
>>
>> I've pretty much got everything working now as it should but there one
>> thing that is throwing me. In the above script k7 was fairly straight
>> forward but suppose I want to calculate the rates for something like:
>>
>> ######################################################
>> import yt
>> from yt.utilities.physical_constants import kboltz
>> import numpy as np
>>
>> def _k17(field, data):
>>     T = data["Temperature"]
>>     print "T.shape = ", T.shape
>>     k17rate = np.zeros(len(T))
>>     i = 0
>>     for t in T:
>>         if t < 1e4:
>>             k17rate[i] = 1.0e-8*pow(t, -0.4)
>>         else:
>>             k17rate[i] = 4.0e-4*pow(t, -1.4)*exp(-15100.0/t)
>>         i = i + 1
>>     return k17rate #cm^3 s^-1
>>
>>
>> yt.add_field("k17rate", units="cm**3/s", function=_k17)
>> Filename = "RD0000/RD0000"
>>
>> ds = yt.load(Filename)
>> sp = ds.sphere('max', (100, 'kpc'))
>> print sp["k17rate"]
>> ######################################################
>>
>> This time the function is messy. There is probably a nice python way of
>> doing the calculation
>> without looping - maybe something similar to k17rate[T < 1e4] =
>> 1.0e-8*np.power(T, -0.4) ?
>>
>> Either way this fails because the array that gets passed to k17 has a
>> shape of (16, 16, 16) and so everything falls
>> apart. The _k17 function actually gets called 5 times in total - firstly
>> with the array of shape (16, 16, 16) and then
>> with flat arrays but it fails on the first non-flat array anyway.
>>
>> Just wondering if you know what's going on here (why the function gets
>> called 5 times and why the first array in a 3D array).
>> Any solutions would of course also be appreciated!
>>
>> Cheers,
>> John
>>
>>
>>
>>
>>
>>
>>
>> On Fri, Sep 26, 2014 at 2:41 PM, John Regan <johnanthonyregan at gmail.com>
>> wrote:
>>
>>> Yip that works! Thanks Britton!
>>>
>>> On Fri, Sep 26, 2014 at 2:35 PM, Britton Smith <brittonsmith at gmail.com>
>>> wrote:
>>>
>>>> Hi John,
>>>>
>>>> I think you should be able to wipe out the units being used in the
>>>> function and add your own by doing:
>>>>
>>>> return data.ds.arr(krate.d, "cm**3/s")
>>>>
>>>> krate.d gives you just the data array without units.
>>>>
>>>> Britton
>>>>
>>>> On Fri, Sep 26, 2014 at 7:05 AM, John Regan <johnanthonyregan at gmail.com
>>>> > wrote:
>>>>
>>>>> Hi All,
>>>>>
>>>>> I'm having a little trouble with Units in YT-3.0.
>>>>>
>>>>> I want to calculate the rates for some collisional dissociations as
>>>>> used in Enzo. The rates are fits so the dimensions are not going to be as
>>>>> required and YT doesn't like that!
>>>>>
>>>>> An example script showing the problem is shown below:
>>>>>
>>>>> #####################################################
>>>>> import yt
>>>>> from yt.utilities.physical_constants import kboltz
>>>>> import numpy as np
>>>>>
>>>>> def _k7(field, data): #H + e -> HM + gamma
>>>>>     Tergs = data["Temperature"]*kboltz
>>>>>     TeV = Tergs.convert_to_units('eV')
>>>>>     krate = 6.77e-15*np.power(TeV, 0.8779)
>>>>>     return krate #cm^3 s^-1
>>>>>
>>>>> yt.add_field(("gas", "k7rate"), units="cm**3/s", function=_k7)
>>>>>
>>>>> Filename = "....../RD0000/RD0000"
>>>>>
>>>>> ds = yt.load(Filename)
>>>>> sp = ds.sphere('max', (100, 'kpc'))
>>>>> print sp["k7rate"]
>>>>> #######################################################
>>>>>
>>>>> The script fails with the error:
>>>>> yt.utilities.exceptions.YTUnitConversionError: Unit dimensionalities
>>>>> do not match. Tried to convert between eV**(8779/10000) (dim
>>>>> (length)**(8779/5000)*(mass)**(8779/10000)/(time)**(8779/5000)) and cm**3/s
>>>>> (dim (length)**3/(time)).
>>>>>
>>>>> which makes sense because the _k7() function is just a fit based on
>>>>> temperature. Is there a way to "set" the units of k7 to be cm^3 s^-1?
>>>>>
>>>>> Cheers,
>>>>> John
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> 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/20140926/8e1d1713/attachment.htm>