[yt-users] Overplotting equipotential contours from total potential
Matthew Turk
matthewturk at gmail.com
Wed Sep 25 04:54:29 PDT 2013
Hi Roberto,
On Tue, Sep 24, 2013 at 9:59 PM, trobolo dinni
<trobolo.trobolo.dinni5 at gmail.com> wrote:
> Hi Matt,
>
> yes, I am using Enzo.
> I am sure that my gravitational field does not include the particles, in
> fact the gravity contribution of their masses is missing from the potential
> (I am not using standard dark matter particles and I think the contribution
> to gravity of this particle type is added after Enzo computes the standard
> potential field). In fact if for example I plot the equipotential contours
> for the default Grav_Potential I get just the equipotential surfaces for the
> gas.
>
> I dug in yt yesterday to understand what was going on and why it did not
> like my potential, in fact I checked exactly the quantities you asked me for
> bot Grav_Potential and my user defined Total_Grav_Potential (in cgs):
>
> In [1309]: grav_pot.min()
> Out[1309]: -38532761511174.227
> In [1310]: grav_pot.max()
> Out[1310]: -1240331838330.0188
> In [1311]: grav_pot.size
> Out[1311]: 16777216
>
> In [1313]: tot_grav_pot.min()
> Out[1313]: -417049888736805.5
> In [1314]: tot_grav_pot.max()
> Out[1314]: -3796359609406.9697
> In [1315]: tot_grav_pot.size
> Out[1315]: 16777216
>
> I checked that the values were changing correctly as the sum of the the
> three partial potentials (gas, part1 and part2), and they are ok. Also the
> total length of the array is ok. So what I did next was defining the
> Total_Grav_Potential field as
>
> def _Total_Grav_Potential(field,data):
>
> #adds the current particle potential to the total potential
> total_grav_potential = data['Grav_Potential']
> return total_grav_potential
>
> so that it just returned the standard Grav_Potential field. Also in this
> case the result was the same, i.e. I got the same error if clim was not
> defined and no contours when clim was defined. Hence the problem was in how
> the field was added to the yt fields.
> I went checking the Creating Derived Fields page in the yt docs and the
> fields.py file from my yt to see how fields were created, so I started
> adding parameters to the add_field() function one at the time to understand
> if they were responsible for the problem.
> Turned out that the validators were the problem, in particular I had to
> validate the following fields:
>
> validators=ValidateDataField(["Grav_Potential","particle_position_x","particle_position_z","particle_position_z"])
>
> which are the exactly the native Enzo fields that I used inside the function
> definition (the used derived fields are instead ok!). So in the end I used
> this:
>
> add_field("Total_Grav_Potential", \
> function=_Total_Grav_Potential, \
> take_log= False , \
> particle_type= True, \
>
> validators=ValidateDataField(["Grav_Potential","particle_position_x","particle_position_z","particle_position_z"]),
> \
> vector_field= False, \
> units=r"\rm{erg}/\rm{g}")
>
> and now I get the contours.
> I uploaded the plots here: http://postimg.org/gallery/5zomtyto/
>
> The problem is solved but I do not really quite understand why the native
> Enzo fields have to be validated.
Hm, me neither. Thank you for your detective work -- this is
something we should look into. Unfortunately today I am otherwise
engaged, but this is something I'd like to look into. Could you file
a bug, target 2.6, and assign to me?
https://bitbucket.org/yt_analysis/yt/issues/new
If you can replicate it with a small script on one of the sample datasets:
yt-project.org/data/
(preferably the IsolatedGalaxy one) that would go a long way toward
helping. The field doesn't even need to mean something, it just needs
to show the behavior you're describing.
> Can I ask you if you can explain me what is happening in your opinion?
>
> Concerning the coding suggestions, looks definitely better to do the
> conversions at the end and use more derived fields, I will rewrite it
> starting from your template.
Let me know if there's anything else we can help with!
-Matt
>
>
> Thank you very much for the help!
> Roberto
>
>
> On 25 September 2013 07:22, Matthew Turk <matthewturk at gmail.com> wrote:
>>
>> Hi Roberto,
>>
>> On Mon, Sep 23, 2013 at 8:59 PM, trobolo dinni
>> <trobolo.trobolo.dinni5 at gmail.com> wrote:
>> > Dear yt users,
>> >
>> > I am running a binary star simulation with enzo, and in my box there is
>> > a
>> > star made up of a gas envelope and a particle as core, and a secondary
>> > star
>> > represented by a single particle (no gas).
>> >
>> > I would like to plot the equipotential surfaces of the system on a slice
>> > plot, hence, since the default Grav_Potential field does not include the
>> > particles potential, I created a custom made field as follows:
>>
>> It looks like you're using Enzo -- is that right? If so, are you sure
>> the Grav_Potential field does not include your particle's potentials?
>> If not, you may be able to apply them as analytic expressions.
>> Regardless, I've supplied some comments below:
>>
>> >
>> > def _Total_Grav_Potential(field,data):
>> > #G in cgs
>> > grav_const_cgs = 6.674e-8
>> > #data for gas
>> > gas_grav_potential_cgs = data['Grav_Potential'] *
>> > (length_unit1/time_unit1)**2.0 #in units of velocity^2
>> > gas_position_x_cgs = data['x'] * length_unit1
>> > gas_position_y_cgs = data['y'] * length_unit1
>> > gas_position_z_cgs = data['z'] * length_unit1
>>
>> You like will want to apply the conversion at the end of your
>> computation. You can also use fields like Radius and ParticleMass to
>> get things back in cm and g, respectively. There are probably better
>> ways than what I am about to propose (Stephen Skory and Elizabeth
>> Tasker have both done things like this) but here's just something off
>> the top of my head.
>>
>> f = data["Grav_Potential"].copy()
>> for i in range(len(data["ParticleMass"])):
>> c = np.array([data["particle_position_%s" % ax] for ax in 'xyz'])
>> data.set_field_parameter("center", c)
>> R = data["Radius"] # an array of gas radii
>> m = data["ParticleMass"][i]
>> # compute grav potential with particle here ...
>> f += ...
>> return f
>>
>> > #initialize the array for the total grav potential to the gas grav
>> > potential
>> > total_grav_potential_cgs = gas_grav_potential_cgs
>> > #loops on the particles
>> > for i in range(len(data['ParticleMassMsun'])):
>> > #data of the current particle
>> > current_particle_mass_cgs= data['ParticleMassMsun'][i] * mass_unit
>> > current_particle_position_x_cgs = data['particle_position_x'][i] *
>> > length_unit1
>> > current_particle_position_y_cgs = data['particle_position_y'][i] *
>> > length_unit1
>> > current_particle_position_z_cgs = data['particle_position_z'][i] *
>> > length_unit1
>> > #computes the array of grav potential for the current particle
>> > current_particle_grav_potential_cgs = -(grav_const_cgs *
>> > current_particle_mass_cgs)/(((gas_position_x_cgs -
>> > current_particle_position_x_cgs)**2.0 + (gas_position_y_cgs -
>> > current_particle_position_y_cgs)**2.0 + (gas_position_z_cgs -
>> > current_particle_position_z_cgs)**2.0)**0.5)
>> > #adds the current particle potential to the total potential
>> > total_grav_potential_cgs = total_grav_potential_cgs +
>> > current_particle_grav_potential_cgs
>> > return total_grav_potential_cgs
>> > #adds the field to the available yt fields
>> > add_field("Total_Grav_Potential", function=_Total_Grav_Potential,
>> > take_log=True ,units=r"\rm{erg}/\rm{g}")
>> >
>> > where the length unit and time unit variables have the function to
>> > convert
>> > my quantities in cgs.
>> > If I then call my "Total_Grav_Potential" field, I get reasonable values,
>> > but
>> > when I try to inspect the shape of the equipotential surfaces on the
>> > xy-plane by plotting the following:
>> >
>> > sp = SlicePlot(pf, 'z', "Density", width = 1.0)
>> > sp.annotate_velocity(factor=16, normalize=True)
>> >
>> > sp.annotate_contour('Total_Grav_Potential',clim=(max_total_pot,min_total_pot))
>> > sp.annotate_particles(1.0, p_size=50.0, marker='o', col='black')
>> > sp.annotate_text((0.7,1.05), "time = "+str(pf.current_time)+"yr")
>> > sp.save(plot_dir+"density_slice_z_"+str(index)+".png")
>> >
>> > nothing happens, and the plot does not show any additional contour;
>> > while if
>> > I do the same with the default "Grav_Potential" field, the gas
>> > equipotential
>> > contours get plotted without problems.
>> > I also tried to play with the "annotate_contour" options "ncont",
>> > "factor"
>> > and "clim" to refine my plot but the result is always the same.
>> > Additionally, if I do not set the "clim" option the following error
>> > shows
>> > up:
>> >
>> > ValueError: zero-size array to minimum.reduce without identity
>>
>> This message suggests something funny may be going on. What is the result
>> of:
>>
>> dd = pf.h.all_data()
>> dd["Total_Grav_Potential"].min(), dd["Total_Grav_Potential"].max(),
>> dd["Total_Grav_Potential"].size
>>
>> >
>> >
>> > I would like to ask if my custom definition of the Total_Grav_Potential
>> > field is correct from a yt coding point of view and if I should modify
>> > it
>> > (or if there is any improvement I can do).
>> > Also, I would like to ask if you have any idea of what could be the
>> > problem
>> > with the plotting.
>>
>> Let us know how it goes, and I hope this helps!
>>
>> -Matt
>>
>> >
>> >
>> > Thanks for the help,
>> > Roberto
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> > _______________________________________________
>> > 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
>
More information about the yt-users
mailing list