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

[Nathan Goldbaum]

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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20160518/366758ff/attachment.html>