[yt-users] Defining a new field, including IF statement

[Britton Smith]

Hi Yusuke,

This error comes from the fact that your field function is being evaluated
for an array of elements all at once.  You will likely have portions of the
array that satisfy both criteria.  You'll need to first define a field
array and then set the values for the portions above the density threshold
and then below.  I recommend doing something like the following:

def _SFR_Density(field, data):
   f_sf = 0.01
   rho_thresh = 57.5 * Mu * hydrogen_mass
   # make field array
   field_data = data["density"].copy()
   # density threshold
   threshold = data["density"] > rho_thresh
   # set values above threshold
   field_data[threshold] *= f_sf * / (3*pi/32/Grav/data["density"])**(0.5)
   # set values below
   field_data[np.invert(threshold] = 0
   return field_data

The threshold array is a boolean array that can be used to access all the
values satisfying your criteria all at once.  Then, the np.invert gives you
all the values not satisfying it.  Hope that helps!

Britton

On Wed, Nov 1, 2017 at 10:46 PM, Yusuke Fujimoto <
yusuke.fujimoto.jp at gmail.com> wrote:

> Dear yt-users,
>
> I am trying to define a new field, SFR_Density, which includes if
> statements.
> That is defined by,
> f_sf * rho / t_ff  if rho > rho_threshold
> 0                      if rho < rho_threshold
> (f_sf is star formation efficiency, rho is gas density, and t_ff is free
> fall time)
>
> Here is my definition I tried.
>
> def _SFR_Density(field, data):
>
>    f_sf = 0.01
>
>    rho_thresh = 57.5 * Mu * hydrogen_mass
>
>    if data["density"] > rho_thresh:
>
>       return f_sf * data["density"] / (3*pi/32/Grav/data["density"])**(0.5
> )
>
>    else:
>
>       return 0.0 * f_sf * data["density"] / (3*pi/32/Grav/data["density"])
> **(0.5)
>
> add_field("SFR_Density", function=_SFR_Density)
>
> But it fails with an error message of "ValueError: The truth value of an
> array with more than one element is ambiguous. Use a.any() or a.all()"
> I am sure this is because data["density"] is a list, so I need to specify
> the exact index number of elements at the if statement. I have tried many
> tests, but now I am in stuck.
>
> Any help would be very appreciated.
>
> Thanks,
> Yusuke
>
> _______________________________________________
> 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/20171102/7f10e86c/attachment-0002.html>