[yt-users] Problem with gradient of a derived field

[Nathan Goldbaum]

Hi José,

You're right, of course. This incorrect behavior is due to a bug in yt.
Thank you for the report! Sorry for being doubtful initially, it took me a
while to convince myself that yt's answers were wrong.

This is a particularly pernicious error, since it will only show up for
fields that need field parameters and ghost zones, which unfortunately is
apparently not well tested...

I've made a pull request to fix this here:

https://bitbucket.org/yt_analysis/yt/pull-requests/2186

With this fix, I now get this plot for the gradient of your Omega field:

http://i.imgur.com/DEo4mBw.png

Let me know if you need help merging in my fix locally. Hopefully it will
be merged into the development branch of yt within the next few days.

-Nathan

On Wed, May 18, 2016 at 9:28 PM, José Utreras <jutreras at das.uchile.cl>
wrote:

> Thanks Nathan! Now I don't see that division in the field.
> Still, I think that it does not have the correct values. In this case,
> with the correct definition of the field, the derivative should be positive
> for x>0 and negative for x<0.
>
> 2016-05-18 20:43 GMT-03:00 Nathan Goldbaum <nathan12343 at gmail.com>:
>
>> Hi Jose,
>>
>> On Wed, May 18, 2016 at 1:17 PM, José Utreras <jutreras at das.uchile.cl>
>> wrote:
>>
>>> Hi everyone,
>>> I'm having the following problem: I want to calculate the gradient field
>>> of the angular velocity which I have defined as follows:
>>>
>>> >def _Omega(field,data):
>>> >        if data.has_field_parameter("bulk_velocity"):
>>> >                bv =
>>> data.get_field_parameter("bulk_velocity").in_units("cm/s")
>>> >        else:
>>> >                bv = data.ds.arr(np.zeros(3), "cm/s")
>>> >        xv = data["gas","velocity_x"] - bv[0]
>>> >        yv = data["gas","velocity_y"] - bv[1]
>>> >        center = data.get_field_parameter('center')
>>> >        x_hat = data["x"] - center[0]
>>> >        y_hat = data["y"] - center[1]
>>> >        r = np.sqrt(x_hat*x_hat+y_hat*y_hat)
>>> >        x_hat /= r**2
>>> >        y_hat /= r**2
>>> >
>>> >        return np.sqrt((xv*y_hat-yv*x_hat)**2)
>>>
>>
>> The issue is this line ^
>>
>> You're taking the absolute value of z component of the angular velocity,
>> which introduces a discontinuity in the first derivative - that's the
>> artifact you're seeing in your plots.  I think you really want to have:
>>
>> > return (xv*y_hat-yv*x_hat)
>>
>> Here's what Omega and the gradient look like with that definition:
>>
>> http://i.imgur.com/7ElXqw6.png
>> http://i.imgur.com/fjoLXX3.png
>>
>>
>>> >yt.add_field("Omega", function=_Omega,take_log=False,
>>> units=r"1/yr",validators=[ValidateParameter('bulk_velocity')])
>>>
>>> This is the profile I have got using the SlicePlot tool:
>>>
>>>
>>> https://drive.google.com/file/d/0B43CkTTj8OeOcy16TkFaZlBhRE0/view?usp=sharing
>>>
>>> So, it is a function of radius. The density is also a function of
>>> radius. The problem is that if I use
>>>
>>> >grad_fields = ds.add_gradient_fields(("gas","Omega"))
>>> >slice = yt.SlicePlot(ds, 'z', "Omega_gradient_x")
>>>
>>> I get this
>>>
>>>
>>> https://drive.google.com/file/d/0B43CkTTj8OeOLUxGVUJSQ253TE0/view?usp=sharing
>>>
>>> ​the sign should only change across x=0. I did the same for the density
>>> field getting the correct behaviour.
>>> I've also tried to define the gradient field
>>>
>>> >def _vorticity_2(field, data):
>>> >    f  = (data[ftype, "Omega"][sl_right,sl_center,sl_center] -
>>> >          data[ftype, "Omega"][sl_left,sl_center,sl_center]) \
>>> >          / (div_fac*just_one(data["index", "dx"]))
>>> >
>>> >    new_field = data.ds.arr(np.zeros_like(data[ftype, "velocity_z"],
>>> >                                          dtype=np.float64),
>>> >                            f.units)
>>> >    new_field[sl_center, sl_center, sl_center] = f
>>> >    return new_field
>>>
>>> But I got the same results.
>>>
>>> If I'm doing something wrong please let me know.
>>>
>>> Cheers!
>>>
>>> Jose Utreras
>>>
>>> _______________________________________________
>>> 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/20160518/e802a32f/attachment.html>