Mailman 3 June 2010 - yt-users

total number of cells
by Stephen Skory 05 Jun '10

05 Jun '10

Hi All, what's the easiest way to get the number of cells in a simulation? Thanks! _______________________________________________________ sskory(a)physics.ucsd.edu o__ Stephen Skory http://physics.ucsd.edu/~sskory/ _.>/ _Graduate Student ________________________________(_)_\(_)_______________

2 1

error Zero sized dimension for non-unlimited dimension
by Jean-Claude Passy 05 Jun '10

05 Jun '10

Dear all, I am encountering some issues when running the script enclosed: ###################################################################################################### loading DM_part/DD0001/CommonEnvelope0001 yt INFO 2010-06-03 17:36:57,270 Storing the binary hierarchy yt INFO 2010-06-03 17:36:57,523 Gathering a field list (this may take a moment.) Traceback (most recent call last): File "script2body.py", line 39, in <module> region = data.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]) File "/rpod2/jcpassy/yt-x86_64/src/yt-trunk-svn/yt/lagos/OutputTypes.py", line 150, in _get_hierarchy self.__hierarchy = self._hierarchy_class(self, data_style=self.data_style) File "/rpod2/jcpassy/yt-x86_64/src/yt-trunk-svn/yt/lagos/HierarchyType.py", line 374, in __init__ AMRHierarchy.__init__(self, pf, data_style) File "/rpod2/jcpassy/yt-x86_64/src/yt-trunk-svn/yt/lagos/HierarchyType.py", line 64, in __init__ self._detect_fields() File "/rpod2/jcpassy/yt-x86_64/src/yt-trunk-svn/yt/lagos/HierarchyType.py", line 615, in _detect_fields self.save_data(list(field_list),"/","DataFields",passthrough=True) File "/rpod2/jcpassy/yt-x86_64/src/yt-trunk-svn/yt/lagos/HierarchyType.py", line 192, in _save_data arr = myGroup.create_dataset(name,data=array) File "/rpod2/jcpassy/yt-x86_64/lib/python2.6/site-packages/h5py-1.2.0-py2.6-linux-x86_64.egg/h5py/highlevel.py", line 353, in create_dataset return Dataset(self, name, *args, **kwds) File "/rpod2/jcpassy/yt-x86_64/lib/python2.6/site-packages/h5py-1.2.0-py2.6-linux-x86_64.egg/h5py/highlevel.py", line 856, in __init__ space_id = h5s.create_simple(shape, maxshape) File "/rpod2/jcpassy/yt-x86_64/lib/python2.6/site-packages/h5py-1.2.0-py2.6-linux-x86_64.egg/h5py/_sync.py", line 88, in wrap return func(*args, **kwds) File "h5s.pyx", line 106, in h5py.h5s.create_simple (h5py/h5s.c:973) h5py._stub.ValueError: Zero sized dimension for non-unlimited dimension (Function arguments: Bad value) ###################################################################################################### The weird thing is that if I run the script again, it works and then breaks when trying to read the next file (DM_part/DD0001/CommonEnvelope0001), and so on... Data contain only 2 particles without gas and as far as I remember, I could execute those commands without problems before. I would appreciate if someone could help me with that. Thanks a lot, Jean-Claude #!/rpod2/jcpassy/yt-x86_64/bin python2.6 import sys import commands as C import numpy as N import pylab as P import yt.raven as R import yt.lagos as L filebase = sys.argv[-1] minSnap = 0 maxSnap = 401 skip = 1 t = N.zeros(maxSnap/skip) X = N.zeros((maxSnap/skip,3,2)) V = N.zeros((maxSnap/skip,3,2)) r = N.zeros(maxSnap/skip) X2 = N.zeros((maxSnap/skip,3,2)) V2 = N.zeros((maxSnap/skip,3,2)) r2 = N.zeros(maxSnap/skip) j = 0 for i in range(minSnap,maxSnap,skip): # File 1 filen = 'DM_part/'+ 'DD' + '%04i' % i + '/' + filebase + '%04i' % i print "loading %s" % filen data = L.EnzoStaticOutput(filen,data_style="enzo_packed_3d") # if no gas t[j] = data["InitialTime"] * data["years"] region = data.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]) i1,i2 = region['particle_index'] X[j,0,i1],X[j,0,i2] = region['particle_position_x'] X[j,1,i1],X[j,1,i2] = region['particle_position_y'] X[j,2,i1],X[j,2,i2] = region['particle_position_z'] V[j,0,i1],V[j,0,i2] = region['particle_velocity_x'] V[j,1,i1],V[j,1,i2] = region['particle_velocity_y'] V[j,2,i1],V[j,2,i2] = region['particle_velocity_z'] r[j] = pow((X[j,0,i1] - X[j,0,i2]),2) + pow((X[j,1,i1] - X[j,1,i2]),2) + pow((X[j,2,i1] - X[j,2,i2]),2) r[j] = pow(r[j],0.5) j = j + 1

2 2

Re: [yt-users] fields in ghost zones
by David Collins 03 Jun '10

03 Jun '10

Hi, Jean-Claude-- (I copied this to the list, so all may enjoy my wisdom.) HDF5 brings some extra work over straight ascii files, but it's 1000% worth it in the long run. It's your call if you think it's right for this problem. I'm the kind of person who will go far out of my way to be lazy later, so this is what I'd do. (Even if it's not right for this problem, HDF5 is worth knowing. I hate all datasets that aren't in HDF5, now that I've gotten comfortable with it.) The HDF5 website is here: http://www.hdfgroup.org/HDF5/ I've been told that the HDF5 Lite interface is pretty good: http://www.hdfgroup.org/HDF5/Tutor/h5lite.html though I haven't used it myself, the syntax is supposed to be simpler. It's worth going through a few examples on the website to get used to it. I would do it using a similar syntax that's in Grid_Group_ReadGrid.C, though this file has a bunch of enzo baggage in it, so it's kind of difficult to read. One thing to note is that HDF5 is a C library, but Enzo arrays are stored in Fortran ordering. So be careful that the array will have the X and Z indices swapped if you go to 3d. Another thing is to watch out for data types. The short version: HDF5 files have a directory structure like a unix directory tree. Enzo output files have a format like File/Grid Group/ Dataset Each processor does it's own output, so there should be one cpu???? file for each MPI task. So you need to do the following in your code: 1.) Open the file (data0000.cpu0001) 2.) Open the "Group", in this case the grid ID, which will be something like Grid00001. h5ls will show you the contents of the file 3.) Open the dataset ('Density', 'x-velocity', etc) 4.) Read the dataset 5.) Close the dataset, group, file Syntax and details can be found on the hdf5 websites or in the Enzo Grid_Group_ReadGrid.C files. Hope that helps! d. On Wed, Jun 2, 2010 at 2:06 PM, Jean-Claude Passy <jcpassy(a)gmail.com> wrote: > Hi David, > > it looks like a good idea, thanks. > Could you just be more specific and explain briefly how to read the HDF5 > file, especially in a C++ program ? I am sorry but I am also not familiar > with the h5ls command. > > Thanks a lot, > > Jean-Claude > > David Collins a écrit : > > Hi, JC-- > > > Would it be easier to skip YT all together, and just have your problem > generator read the HDF5 file produced by the 1d run directly? This would > streamline your processes, and you wouldn't have the noise incurred by both > the > ASCII file format and the two units conversions. > > I would probably approach your task by looping over only the active > zones in enzo, and incrementing an extra counter. > > int counter = 0; > for(int field=0;field<NumberOfBaryonFields;field++) > for(int k=GridStartIndex[2];k<=GridEndIndex[2]; k++) > for(int j=GridStartIndex[1];j<=GridEndIndex[1]; j++) > for(int i=GridStartIndex[0];i<=GridEndIndex[0]; i++){ > index = i + GridDimension[0]*(j + GridDimension[1]*k); > BaryonField[ DensNum ][ i ] = DensityFromFile[ counter ]; > counter++ > } > > > > This code snippit came from the wiki page on Baryon access, for future > reference: > http://lca.ucsd.edu/projects/enzo/wiki/Tutorials/BaryonFieldAccess > > Does that help? > d. > > > On Wed, Jun 2, 2010 at 12:07 PM, Jean-Claude Passy <jcpassy(a)gmail.com> > wrote: > > > Hi Matt, > > I will try to explain what my goal is. > I have a 1D profile of a star that is physically at the equilibrium. > Regarding the problem initialization, I set up the density and the total > specific energy fields according to the model, velocities equal 0. However, > the star is not numerically stable so I have to relax it. Therefore, I > divide the velocity by 2 at each timestep, and I let Enzo run for a few > dynamical times. Then, I read the last data dump and create an > relaxed_model.dat: > > ########################################################## > region = pf.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]) > output = 'relaxed_model.dat' f = open(output,'w') print >>f, '# Grid Index > rho pressure' density = region["Density"] pressure = region["Pressure"] size > = NP.size(density) > for i in range(0,size,1): dtmp = density[i] ptmp = > pressure[i] print >>f, "0", i, '%06e' % dtmp, '%06e' % ptmp > ########################################################## > > so I end up with a file like this: > > ########################################################## # Restart file > for only one grid > # Grid Index rho pressure 0 0 5.099011e-10 2.568020e+05 0 1 5.099011e-10 > 2.568020e+05 0 2 5.099011e-10 2.568020e+05 > ... > ########################################################## > Finally, I read index, density, pressure from that previous file and do: > > ########################################################## for (k = 0; k < > GridDimension[2]; k++) for (j = 0; j < GridDimension[1]; j++) > for (i = 0; i < GridDimension[0]; i++) { > index = i + GridDimension[0]*(j + GridDimension[1]*k); > density /= DensityUnits; pressure /= PressureUnits; > BaryonField[0][index] = density; BaryonField[1][index] = pressure > / ((Gamma - 1.0) * density); > ########################################################## > The problem is GridDimensions contain the ghost zones so the variable index > does not match with the actual index read in the file relaxed_model.dat. > That is why I wanted to have that file values for the ghost zones as well. > > Does it make sense ? Do you have any suggestion ? > > Thanks for your help, > > > JC > > Matthew Turk a écrit : > > Hi Jean-Claude, > > There are a couple aspects to this. The first is that Enzo doesn't > output the ghost zones -- so any ghost zones handled inside yt are > generated by yt. Were Enzo to output the ghost zones, we would > probably be able to handle this, but it does not. > > Derived fields can depend on the generation of ghost zones, but keep > in mind that these are ghost zones generated by yt. These ghost zones > are constructed in a similar manner to how Enzo generates them, but > there may be minor differences. You can manually inspect ghost zones > on a *grid* by calling retrieve_ghost_zones on that grid. > > If you could tell us a bit more about your goal, maybe we could help > out a bit more? > > -Matt > > On Wed, Jun 2, 2010 at 9:32 AM, Jean-Claude Passy <jcpassy(a)gmail.com> wrote: > > > Dear all, > > in order to set up my simulation, I need to access the values of certain > fields in the ghost zones. > If I do something like: > > pf = load(data) region = pf.h.region(...) x = region[field] > > x contains the values of field for the physical grid only. Is there a way I > can get the same 1D-array but with the ghost zones included as well ? > > Thanks for your help, > > Jean-Claude > > _______________________________________________ > yt-users mailing list > yt-users(a)lists.spacepope.org > http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org > > > > > _______________________________________________ > yt-users mailing list > yt-users(a)lists.spacepope.org > http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org > > > > > -- Sent from my Stone Tablet and carried by my Pterodactyl.

2 1

Profiles and zero-valued fields
by Matthew Turk 03 Jun '10

03 Jun '10

Hi all, After a chat with David Collins today, I wanted to share a relatively quick note about the profiles and how they behave with fields, specifically the field "Convergence", but any fields that zero out values are equally applicable -- this will also extend to deposited particle fields. The issue here is in how the profiling works. The 1- and 2-D profile objects are computed as: sum( f_i * w_i ) / sum( w_i ) Note that for fields where the inapplicable values are zero'd out -- in this case, "Convergence" zeroes out all positive divergence calculations and returns only the negative values, but sign-flipped. However, the same is not done for the denominator! So your profile will give you the "what you said" results but not the "what you mean" results. The solution here is to use either an extracted region or an inline extracted region (for large datasets): dd = pf.h.all_data() convergent_region = dd.extract_region( dd["Divergence"] < 0.0 ) pc.add_profile_object(convergent_region, ["Density", "Convergence"]) -or- dd = pf.h.all_data() convergent_region = dd.cut_region( [ "grid['Divergence'] < 0.0" ]) pc.add_profile_object(convergent_region, ["Density", "Convergence"]) Anyway, I hope that if anybody else runs into something similar that this will help out! Best, Matt

1 0

issue with wx python and matplotlib
by Eric Hallman 03 Jun '10

03 Jun '10

Hello, I am having a problem with yt on ubuntu 9.10 that I haven't seen before, wondering if anyone else has this error. I have wxPython 2.8.10 installed, so the error message makes no sense to me. Thanks for any help. Python 2.6.3 (r263:75183, Jan 29 2010, 11:07:13) [GCC 4.4.1] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> from matplotlib.pyplot import * Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/home/hallman/yt/yt-unknown/lib/python2.6/site-packages/ matplotlib/pyplot.py", line 78, in <module> new_figure_manager, draw_if_interactive, show = pylab_setup() File "/home/hallman/yt/yt-unknown/lib/python2.6/site-packages/ matplotlib/backends/__init__.py", line 25, in pylab_setup globals(),locals(),[backend_name]) File "/home/hallman/yt/yt-unknown/lib/python2.6/site-packages/ matplotlib/backends/backend_wxagg.py", line 23, in <module> import backend_wx # already uses wxversion.ensureMinimal('2.8') File "/home/hallman/yt/yt-unknown/lib/python2.6/site-packages/ matplotlib/backends/backend_wx.py", line 137, in <module> raise ImportError(missingwx) ImportError: Matplotlib backend_wx and backend_wxagg require wxPython >=2.8 >>> Eric Hallman Google Voice: (774) 469-0278 hallman13(a)gmail.com

2 1

Cloud in cell deposition
by Cameron Hummels 02 Jun '10

02 Jun '10

Hello everyone: I would like to make some images of the stellar (or DM) particles for some of my enzo halos. YT's current abilities for overplotting particle data on the grid gas data leaves something to be desired (although a long-term solution is being coordinated), so I was looking for an alternative. Matt mentioned that some of you have circumvented this very problem by depositing particles into nearby grid cells and then using YT's excellent rendering abilities to produce some knockout images. Does anyone have any suggestions on how one might go about generating such an image for, say, star particles? Thanks! Cameron -- Cameron Hummels PhD Candidate, Astronomy Department of Columbia University Public Outreach Director, Astronomy Department of Columbia University NASA IYA New York State Student Ambassador http://outreach.astro.columbia.edu

2 1

field gradients / vorticity / ghost cells
by Andrew Cunningham 02 Jun '10

02 Jun '10

I'd like to make a derived field that is the gradient of a primitive field. Is there an example of how I can do this? I'm not sure how to handle ghost cells for something like this. I spoke with Jeff and he indicated that there was a vorticity example somewhere, but I cannot find it. Thanks.

6 18