Overplotting equipotential contours from total potential
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:
* 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* * * * #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* * *
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.
Thanks for the help, Roberto
Hi Roberto,
On Mon, Sep 23, 2013 at 8:59 PM, trobolo dinni trobolo.trobolo.dinni5@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@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
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. 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.
Thank you very much for the help! Roberto
On 25 September 2013 07:22, Matthew Turk matthewturk@gmail.com wrote:
Hi Roberto,
On Mon, Sep 23, 2013 at 8:59 PM, trobolo dinni trobolo.trobolo.dinni5@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@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Hi Roberto,
On Tue, Sep 24, 2013 at 9:59 PM, trobolo dinni trobolo.trobolo.dinni5@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@gmail.com wrote:
Hi Roberto,
On Mon, Sep 23, 2013 at 8:59 PM, trobolo dinni trobolo.trobolo.dinni5@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@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
participants (2)
-
Matthew Turk -
trobolo dinni