[yt-users] Overplotting equipotential contours from total potential

Matthew Turk matthewturk at gmail.com
Tue Sep 24 14:22:29 PDT 2013


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
>



More information about the yt-users mailing list