Mailman 3 October 2013 - yt-users

Flux through spherical surfaces
by Morgan MacLeod Oct. 28, 2013

Oct. 28, 2013

Dear all, I am hoping to calculate the flux of mass ("mdot") and angular momentum ("Ldot") through a series of concentric spherical surfaces in a simulation snapshot (centered around an absorbing "sink" within the grid). So far, I've attempted this along these lines: pf = load(fn) # load data rad = 0.1 sp = pf.h.sphere([0,0,0],1.1*rad) surf = pf.h.surface(sp,"Radius",rad) flux = surf.calculate_flux("x-velocity","y-velocity","z-velocity","Density") where flux gets assigned to the mass accretion rate in this case. I believe I could do the angular momentum by calculating the flux of a derived variable that is defined as data["SpecificAngularMomentumZ"] * data["Density"] such that the resulting units are angular momentum per volume. I was wondering if there is a better way to go about this calculation, it seems that others must have tried to look at similar questions. In particular, defining an "isoradius" surface seems to accomplish what I was hoping to do, but there must be a cleaner way. Thanks so much for your advice, Morgan

2 2

How to use if, elif, else statements in creating a derived field
by Reju Sam John Oct. 27, 2013

Oct. 27, 2013

Dear all, I would like to create a derived field which should return a value according to specified conditions. My definition of new field is shown below.. But it is giving error. @derived_field(name = "MyField") def my_new_field(field, data): if data.pf["MachNumber"] < 2 : return (.5* (data["MachNumber"])**2) - 1 elif data.pf["MachNumber"] < 100 : return (.5* (data["MachNumber"])**4) + (.5* (data["MachNumber"])**3) else: return 5.446361E-01 Please suggest me how to implement if, elif, else statements in creating a derived field. Thank You Reju Sam John

6 12

running out of memory in a loop
by Stephanie Tonnesen Oct. 23, 2013

Oct. 23, 2013

Hi everyone, Okay, so I am trying to make a SlicePlot, a ProjectionPlot, and a phase_object for a bunch of different outputs. My code keeps using up all my memory and dying after the first loop. I am not using variables, I think, so how do I free up my memory at the end of each loop? See code below. Thanks, Stephanie import matplotlib.pylab as pl from mpl_toolkits.mplot3d import Axes3D from yt.mods import * def _xkpc(field,data): return(data["x"]*26.) def _ykpc(field,data): return(data["y"]*26.) def _zkpc(field,data): return(data["z"]*26.) def _CellMassSolar(field,data): return(data["CellMass"]/1.99e33) add_field("xkpc",function=_xkpc,units="kpc") add_field("ykpc",function=_ykpc,units="kpc") add_field("zkpc",function=_zkpc,units="kpc") add_field("CellMassSolar",function=_CellMassSolar,units=r'M$_{\odot}$') outstring = ['0020','0030','0040','0050','0060','0070','0080','0090','0100','0110','0120','0130'] i = 0 while i < len(outstring): pf = load("id0/rps."+outstring[i]+".vtk",parameters={"TimeUnits":3.086e14,"LengthUnits":8.0236e22,"DensityUnits":9.999e-30,"VelocityUnits":2.6e8}) p = SlicePlot(pf,"x","Density") p.set_log("Density", True) p.save() p2 = ProjectionPlot(pf,"x","Density") p2.set_log("Density", True) p2.save() alld = pf.h.all_data() pc = PlotCollection(pf) tracerp25 = alld.cut_region(["grid['specific_scalar[0]'] > 0.25", "grid['z'] > 0.04"]) tp25extrema = tracerp25.quantities["Extrema"]("Density") pc.add_phase_object(tracerp25,["zkpc","z-velocity","CellMassSolar"],weight=None,x_bins=50,y_bins=50,y_bounds = [-1e8,1.5e8],x_log=False,y_log=False) pc.set_zlim(10000,1e10) pc.save() i = i+1

4 6

Changing the cell range in add_phase_object
by Stephanie Tonnesen Oct. 23, 2013

Oct. 23, 2013

Hi everyone, So I am making a phase_object, plotting mass as a function of z vs z_velocity. Some of my [z,z-velocity] bins have very little mass in them, and I would like to limit the range of CellMass that is plotted--right now most of the interesting details are lost in a relatively small range of CellMass. Right now my range in the colormap is 1e35 to 1e42, and I want it to be 1e37 to 1e42, and ignore all bins below 1e37. Can I do that? Thanks, Stephanie

2 3

Re: [yt-users] plotting average quantaties
by Alexander Philippov Oct. 23, 2013

Oct. 23, 2013

Thanks all for very constructive suggestions! Matt and Jeff thing works good. Thanks, Sasha On Tue, Oct 22, 2013 at 5:45 AM, j s oishi <jsoishi(a)gmail.com> wrote: > I think Sasha may want a simpler, unweighted average. You can do this by > setting the weight field to Ones. If you set the weight to None, you'll get > a sum across x. > > You'll probably want to set n_bins to the number of x cells to get what > you are looking for. > > Hope this is helpful, > > Jeff > On Oct 22, 2013 8:25 AM, "Matthew Turk" <matthewturk(a)gmail.com> wrote: > >> Hi Sasha, >> >> I'm not sure I understand what you're looking for. The functionality to >> take averages of one field wrt another field might work for what you're >> looking for, using the BinnedProfile1D object. If so, you could bin by the >> "x" field and then take a weighted average of Density. This might look >> something like this: >> >> dd = pf.h.all_data() >> prof = BinnedProfile1D(dd, n_bins, "x", pf.domain_left_edge[0], >> pf.domain_right_edge[0], log_space = False) >> prof.add_fields(["Density"], weight="CellMassMsun") >> >> Does that help? >> >> Best, >> >> Matt

1 0

plotting average quantaties
by Alexander Philippov Oct. 23, 2013

Oct. 23, 2013

Hi all, I would like to plot 1D profile of quantity (say, density) averaged over 1 axis (say, x) from two-dimensional data. Is there any cheap way to do it? Thanks in advance, Sasha

5 4

Moving a volume rendering camera tutorial
by Alex Bogert Oct. 22, 2013

Oct. 22, 2013

Hello, I'm getting some errors when trying this tutorial : http://yt-project.org/doc/cookbook/complex_plots.html#moving-a-volume-rende… The error: http://pastebin.com/m17u7iCg Thanks, Alex

2 1

Forcing a linear axis (well, actually two) on a phase plot
by Stephanie Tonnesen Oct. 22, 2013

Oct. 22, 2013

Hi Everyone, I am trying to make a phase plot, and I want both of my axes to be linear. I write: pc = PlotCollection(pf) tracerp25 = alld.cut_region(["grid['specific_scalar[0]'] > 0.25", "grid['z'] > 0.01"]) pc.add_phase_object(tracerp25,["z","velocity_z","CellMass"],weight=None) p.set_log("z",False) p.set_log("velocity_z",False) pc.save() And get the figure below. I am betting there is a simple command I am missing. Thanks! Stephanie

2 2

Stars particles
by Francia Riesco Oct. 22, 2013

Oct. 22, 2013

Hello Everybody, I hope you are great! I have a question for you I am analyzing an Enzo data trying to get the center of mass of the Stars Particles. I am suing FOF linking_length = 0.017094633 haloes = FOFHaloFinder(pf, link=linking_length ,dm_only=False) for halo in haloes: ct = halo["creation_time"] sm = halo["ParticleMassMsun"] metal = halo["metallicity_fraction"] stars = (ct > 0) ct = ct[stars] sm = sm[stars] Is there a way to get sm.center_of_mass() or something similar to get the Center of Mass of the Stars Particles?? Thank you in Advance !

2 1

Athena thing Re: yt-users Digest, Vol 68, Issue 29
by Alexander Philippov Oct. 21, 2013

Oct. 21, 2013

Hi John, looks like my version of join vtk script was broken. Stefaphanie Tonnesen gave me a working one, so now it works good. Sasha

1 0