<div dir="ltr">Dear yt users,<div><br></div><div>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).</div>
<div><br></div><div>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:</div>
<div><br></div><div><i> <span class="" style="white-space:pre"> </span>def _Total_Grav_Potential(field,data): </i></div><div><i><span class="" style="white-space:pre"> </span>#G in cgs</i></div><div><i><span class="" style="white-space:pre"> </span>grav_const_cgs = 6.674e-8</i></div>
<div><span class="" style="white-space:pre"><i> </i></span></div><div><i><span class="" style="white-space:pre"> </span>#data for gas</i></div><div><i><span class="" style="white-space:pre"> </span>gas_grav_potential_cgs = data['Grav_Potential'] * (length_unit1/time_unit1)**2.0 #in units of velocity^2</i></div>
<div><i><span class="" style="white-space:pre"> </span>gas_position_x_cgs = data['x'] * length_unit1</i></div><div><i><span class="" style="white-space:pre"> </span>gas_position_y_cgs = data['y'] * length_unit1</i></div>
<div><i><span class="" style="white-space:pre"> </span>gas_position_z_cgs = data['z'] * length_unit1</i></div><div><span class="" style="white-space:pre"><i> </i></span></div><div><i><span class="" style="white-space:pre"> </span>#initialize the array for the total grav potential to the gas grav potential</i></div>
<div><i><span class="" style="white-space:pre"> </span>total_grav_potential_cgs = gas_grav_potential_cgs</i></div><div><span class="" style="white-space:pre"><i> </i></span></div><div><i><span class="" style="white-space:pre"> </span>#loops on the particles</i></div>
<div><i><span class="" style="white-space:pre"> </span>for i in range(len(data['ParticleMassMsun'])):</i></div><div><span class="" style="white-space:pre"><i> </i></span></div><div><i><span class="" style="white-space:pre"> </span>#data of the current particle</i></div>
<div><i><span class="" style="white-space:pre"> </span>current_particle_mass_cgs= data['ParticleMassMsun'][i] * mass_unit</i></div><div><i><span class="" style="white-space:pre"> </span>current_particle_position_x_cgs = data['particle_position_x'][i] * length_unit1</i></div>
<div><i><span class="" style="white-space:pre"> </span>current_particle_position_y_cgs = data['particle_position_y'][i] * length_unit1</i></div><div><i><span class="" style="white-space:pre"> </span>current_particle_position_z_cgs = data['particle_position_z'][i] * length_unit1</i></div>
<div><span class="" style="white-space:pre"><i> </i></span></div><div><i><span class="" style="white-space:pre"> </span>#computes the array of grav potential for the current particle</i></div><div><i><span class="" style="white-space:pre"> </span>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)</i></div>
<div><span class="" style="white-space:pre"><i> </i></span></div><div><i><span class="" style="white-space:pre"> </span>#adds the current particle potential to the total potential</i></div><div><i><span class="" style="white-space:pre"> </span>total_grav_potential_cgs = total_grav_potential_cgs + current_particle_grav_potential_cgs</i></div>
<div><span class="" style="white-space:pre"><i> </i></span></div><div><i><span class="" style="white-space:pre"> </span>return total_grav_potential_cgs</i></div><div><i><span class="" style="white-space:pre"> </span>#adds the field to the available yt fields</i></div>
<div><i><span class="" style="white-space:pre"> </span>add_field("Total_Grav_Potential", function=_Total_Grav_Potential, take_log=True ,units=r"\rm{erg}/\rm{g}")</i></div><div><i><br></i></div><div>where the length unit and time unit variables have the function to convert my quantities in cgs.</div>
<div>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:</div><div><br></div><div>
<div><i><span class="" style="white-space:pre"> </span>sp = SlicePlot(pf, 'z', "Density", width = 1.0)</i></div><div><i><span class="" style="white-space:pre"> </span>sp.annotate_velocity(factor=16, normalize=True)</i></div>
<div><i><span class="" style="white-space:pre"> </span>sp.annotate_contour('Total_Grav_Potential',clim=(max_total_pot,min_total_pot))</i></div><div><i><span class="" style="white-space:pre"> </span>sp.annotate_particles(1.0, p_size=50.0, marker='o', col='black')</i></div>
<div><i><span class="" style="white-space:pre"> </span>sp.annotate_text((0.7,1.05), "time = "+str(pf.current_time)+"yr") </i></div><div><i><span class="" style="white-space:pre"> </span>sp.save(plot_dir+"density_slice_z_"+str(index)+".png")</i></div>
</div><div><br></div><div>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. </div>
<div>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.</div><div>Additionally, if I do not set the "clim" option the following error shows up:</div>
<div><br></div><div><i>ValueError: zero-size array to minimum.reduce without identity</i><br></div><div><i><br></i></div><div><br></div><div>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).<br>
</div><div>Also, I would like to ask if you have any idea of what could be the problem with the plotting.<br></div><div><br></div><div><br></div><div>Thanks for the help,</div><div> Roberto</div>
<div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div></div>