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

[José Utreras]

Hi Nathan,

Thanks a lot! I have modified that file and is working perfectly. No
problem about being doubtful.

Cheers!

2016-05-19 0:20 GMT-03:00 Nathan Goldbaum <nathan12343 at gmail.com>:

> 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
>>
>>
>
> _______________________________________________
> 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/20160519/6c33b2a4/attachment.html>