Mailman 3 February 2012 - yt-dev

RFC: (yt-3.0) Accessing fluids and particles
by Matthew Turk Feb. 27, 2012

Feb. 27, 2012

Hi all, During the geometry handling hangout today (which I think was quite successful!) the question came up of how to handle accessing particles and fluids came up. The suggestion has been put forth to change all access of field values into .fluid/.field and .particles, with the optional type specifier as the first attribute. This would look like: dd = pf.h.all_data() dd.fluid["Dust"]["Density"] or dd = pf.h.all_data() dd.particles["dark_matter"]["velocity-x"] I have an opinion on what we should do, but I'd like to hear from others before continuing down this path. Here are the options: 1) Keep everything as is for now, requiring users to manually selection things like particle type and fluid type (in case of multiple fluids) and just say who cares if people have multiple particle types. 2) Move to requiring dict-of-dicts for particles only, and allow fluids to exist as is. i.e., dd["Density"] dd.particles["dark_matter"]["velocity-x"] 3) Move to requiring both .fluids and .particles, and make both dict-of-dicts. dd.fluid[:]["Density"] # for all densities summed dd.particles[:]["velocity-x"] 4) Leave things mostly as is, but use tuple-access for types. This would mean: dd["Density"] # gives the total density dd["dust","Density"] # only dust density dd["particle-velocity-x"] dd["PopII", "particle-velocity-x"] 5) Use the object for fluids with optional type specifier, and dict-of-dicts for particles: dd["Density"] dd["dust","Density"] dd.particles["dark_matter"]["velocity-x"] -- Votes? To be honest, thinking about this worries me a bit. I'm inclined not to break old code unless absolutely necessary. But, I would like to support two-fluids in a more reasonable way, and I also want to support multiple particle types (a more common use case, I think) more elegantly. As a rule, even though this is coming down the road a bit, I'd like to try to keep in mind that big changes are frustrating for people, so for each change that is made a good reason for making it should be present. In this case, I initially wanted to make it much easier to select individual particle types, but then the discussion grew from there. We need to consider both aesthetic and technical challenges. I'd like to hear any suggestions or feedback on these items, or alternate suggestions. Thanks, Matt

8 10

RFC: (yt-3.0) Geometry handling redesign
by Matthew Turk Feb. 27, 2012

Feb. 27, 2012

Hi all, This email is going to be long, but the changes I am proposing in it will require feedback. Next Monday at 1PM Eastern I'll be holding a hangout on Google+, if you would like to join in and discuss it with lower latency. (Of course, if nobody shows up, that's okay too.) We don't really have a detailed review process for design decisions like this, and I hope that this discussion will guide developing that kind of review process in the future. Everything described here is in the branch geometry_handling at https://bitbucket.org/MatthewTurk/yt-refactor = Background = Geometry in yt has always been designed around two things: 1) Grid patches 2) Spatial location based on bounding boxes This has presented a number of excellent opportunities for optimizing patch-based AMR data. However, it has posed challenges for supporting particle-based datasets. (Particles inside patch-based data are somewhat easier, although if you want to see one of the worst pieces of yt, look no further than FakeGridForParticles, which mocks up particles as grids with varying x,y,z, so that existing fluid-selection routines can be used.) Furthermore, the state of Octree codes is not good. For relatively small (up to about 64) domains in RAMSES, we can do analysis with relative ease -- volume rendering, projections, fluid analysis, etc, all work very nicely. However, what is required is a costly step of regridding octs at the outset, to form grid patches with some efficiency. This regridding is not stored between runs, and it's not *efficient*. This has worked so far, but it is showing its strain. As it stands, yt cannot directly support SPH codes without somehow making SPH look like grids. Octree codes can't be supported without making them look like grid codes. Both of these are very slow operations. And on top of that, the current state of how data is read in, even for grid codes, is a bit gross. There are a number of conditionals, mechanisms for striding over data (i.e., "_get_grids()" which is called in parallel, and iterates over grids, and sometimes you have to call _get_point_indices and sometimes you don't), and probably worst of all a lot of unnecessary memory duplication because arrays in regions/spheres/etc are constructed by iterating over grids, filling in a list of arrays, and then concatenating that list into a second array. Furthermore, several times recently the question has come up of handling and supporting more generic geometries: spherical, polar, cylindrical. Unfortunately, because we're so heavily based around the idea of grids as is, these are difficult to support. All of our collision detection has been designed with rectilinear surfaces in mind. That being said, I'm really quite happy with how yt *functions* for patch-based grids. I'm just not sure that as it stands it is maintainable, and I think that addressing the geometry issues could ultimately also help out a great deal with scalability. (More on that below.) = Geometry Overhaul = For the last couple weeks, I have been working on an overhaul of the geometry system in yt, attempting to handle parallel striding in a format-agnostic way, deal with generating fields in as simple a way as possible, and also attempting to speed things up. I'm going to summarize below the changes that I've made, the new design, and the current status. == Design: Selectors == As it stands, yt data objects (such as Region, Sphere, etc) are mostly thought of as conduits for data. They apply a criteria to the domain and select all of the fluid points or particles within a region or that fit a specific non-spatial criteria. I have extended this concept to construct very lightweight objects called selectors, similar to how ParticleIO is now conducted for Enzo. These are objects in Cython that accept either a set of bounding-box arrays or a set of spatial points and decide if those points are included in a given region. This is essentially a combination of the get_point_indices and get_cut_mask methods that currently exist. These have all been designed to be GIL-free, so that in the future they can be threaded using OpenMP (although it turns out not very much time is spent on them.) To date I have implemented spheres, regions, disks, cutting planes, slices, ortho rays, and grid collections. They all give bitwise identical results to the old method and are almost all substantially faster. (Cutting planes are particularly better, regions are slightly slower, and the others are all about 25-40% faster.) One primary advantage is that these can (and do) trivially provide counting methods, masking methods, and point-evaluation methods. So we can use the same routines to select particles and individual points (for collision detection) that we use for selecting grid cells. This will allow reuse of data object code with particle outputs, and is a primary advantage of the refactor. This also led to separating selection-based data containers from constructed data containers. == Design: Ditching the Hierarchy == (Hopefully that heading caught some eyes.) In line with attempting to reduce Enzo-isms in the code, I have been moving away from specifically calling something a "hierarchy," and instead moving toward referring to geometry handlers. grids, grid_left_edges, grid_right_edges, and so on will no longer be guaranteed to appear hanging off these objects. The new class structure: GeometryHandler GridGeometryHandler(GeometryHandler) OctreeGeometryHandler(GeometryHandler) and eventually, other methods. The concept behind this is to isolate those items that *must* require domain knowledge, mesh knowledge, or particle/fluid distribution knowledge, and place those in the GeometryHandler. All items specific to patch based or block based AMR will go in GridGeometryHandler and so on. This will allow for fundamentally different method of IO for octrees from grids, for instance. The geometry handlers are all required to implement the following methods: _read_selection(self, fields, data_object, chunk) _chunk(self, data_object, chunking_style, ngz = 0) I will describe chunks later. Note that the first routine, _read_selection, accepts the names of the fields, and the data object that describes how to select which portions of a dataset participate in a given selection. For the specific case of Enzo (as I am still in design phase, I have only been working with Enzo data, although I will move on to FLASH and Boxlib next) read_selection calls a routine that counts and masks grids, batches operations to files, and then returns only the necessary data for a selection. The data object does no discarding of out-of-bounds or child-masked data. == Design: Chunking == Currently, yt uses de facto chunking of data in the form of grids. When Derived Quantities want to handle individual grids, one at a time, they "preload" the data from whatever grids the ParallelAnalysisInterface thinks they deserve. These grids are iterated over, and handled individually, then the result is combined at the end. Profiles do something similar. However, both of these are de facto, and not really designed. They rely on calls to semi-private functions on data objects, manually masking data, on and on. To get around this ugly spot, I have started working on redesigning IO and data handling to function in "chunks." In this method, data objects will request that the geometry handler supply a set of chunks (the function _chunk, described above, conducts this task.) Chunks are of the form (IO_unit, size), where IO_unit is only ever managed or handled by _read_selection, above. To start the Chunks shuffling over the output, the code calls data_source.chunks(fields, chunking_style). Right now only "spatial", "io" and "all" are supported for chunking styles. This corresponds to spatially-oriented division, IO-conserving, and all-at-once (not usually relevant.) What happens next is kind of fun. The chunks function looks like this: def chunks(self, fields, chunking_style, **kwargs): for chunk in self.hierarchy._chunk(self, chunking_style, **kwargs): with self._chunked_read(chunk): self.get_data(fields) yield self Note what it does here -- it actually yields *itself*. However, inside the chunked_read function, what happens is that the attributes corresponding to the size, the current data source, and so on, are set by the geometry handler (still called a hierarchy here.) So, for instance, execution might look like this: for ds in my_obj.chunks(["Density"], "spatial"): print ds is my_obj print ds["Density"].size The first line will actually print True, but the results from the second one will be the size of (for instance) the grid it's currently iterating over. In this way, it becomes much easier to stride over subsets of data. Derived quantities now look like this: chunks = self._data_source.chunks([], chunking_style="io") for ds in parallel_objects(chunks, -1): rv = self.func(ds, *args, **kwargs) It chunks things, evaluates and then stores intermediate results. This is not meant to replace spatial decomposition in parallel jobs, but it *is* designed to enable much easier and *mesh-neutral* division of labor for parallelism and for IO. If we were to call chunk on an octree, it no longer has to make things look like grids; it just makes them look like flattened arrays (unless you chunk over spatial, which I haven't gotten into yet.) Essentially, by making the method of subsetting and striding over subsetted data more compartmentalized, the code becomes more clear and more maintainable. = Status = As it stands, I am able to select data very nicely -- this means slices, profiles, derived quantities all work. However (see below) particles still need some work and handling, and constructed data objects (smoothed covering grids, projections) have yet to be ported. = Open Questions = However, the main reason I am writing -- when the code is not yet fully functional! -- is to ask for advice and consideration on a few points. == Particle Types == In the past, we've talked about handling particle types differently depending on the particle. Because this is destined for yt-3.0, I would like to propose that we break compatibility with old scripts in favor of clarity in our particle types. Specifically, rather than doing: obj["particle_position_x"][obj["particle_type"] == 1] We do, obj[("particle_positions_x", 1)] (or obj[("particle_positions_x", "PopIICluster")] if you are using a code that supports such expression.) And if a field is found to be a particle type field, but without a tuple, we simply raise a KeyError. We can emulate old behavior with: obj[("particle_position_x",None)] But I am inclined to say that I would rather *mandate* particle type selection, rather than merely allow it. This will lead to much greater clarity, I think. [+-][01] on == Particle Access == Currently, we access particles off a data_object: sphere["particle_position_x"] Do we want to keep this, or move to being much more explicit, that we require sphere.particles["particle_position_x"] ? I can see both sides to this. Explicit is good, and would allow us to be both more clear with users and about fields (and with our IO methods) but allowing old behavior encourages not having to necessarily worry about if something is a fluid or a particle, which is a good thing I think. == User Survey == Should we conduct a brief (only a handful of questions) polling users for their feelings on the changes or potential changes to nomenclature, particles, and so on? We could use it also as a time to identify other outstanding issues with yt. = Conclusion = Does this all make a sort of sense? [+-][01] on the design presented so far If you're interested in helping out, please get in touch. I'm still laying groundwork in some areas, but in other areas I would appreciate help and advice moving forward. Thanks to everyone who made it this far, Matt

5 6

.quantities issue
by Stephen Skory Feb. 26, 2012

Feb. 26, 2012

Hi all, when I was trying to work on Mike's latest halo finder problem, I tried swapping out the total_mass sum for a .quantities on this line: https://bitbucket.org/yt_analysis/yt/src/73fa0ace10e3/yt/analysis_modules/h… but doing that seriously changed the answers. Digging a bit deeper I see that for some reason, the same value is being fed in twice (I'm using two cores) to mpi_allreduce, which gives the wrong total_mass out. I see different values going into mpi_allreduce, and the correct total_mass, when I use the original summing method. This has led me to conclude that the .quantities is not being done correctly here - there's something wrong with the multi-level parallelism here I think... or something. Does anyone have any ideas of what could be going on? Thanks! -- Stephen Skory s(a)skory.us http://stephenskory.com/ 510.621.3687 (google voice)

2 5

yt & Enzo w. Zeus
by Elizabeth Tasker Feb. 25, 2012

Feb. 25, 2012

Hi, I wanted to discuss a possibly more consistent way of dealing with velocities in yt when analysing a simulation with Enzo's Zeus hydro scheme (HydroMethod = 2). Because Zeus' velocities are face-centred, not cell centred, there is a half-cell off-set that the routines in yt sometimes allow for and sometimes do not. For instance, the DivV does allow for Zeus' off-set but the calculation of the SpecificAngularMomentum does not. I often hit this problem when examining objects with yt. Either I have to pretend Zeus has cell-centred velocities or write a new routine which is a nuisance. I was wondering if we could introduce a universal field that cell-centred Zeus' velocities. I think it would look like this: http://paste.yt-project.org/show/2192/ But, I have three questions: (1) Does the scheme (using the slice(None, -1, None) for arrays) work if "data" is not a continuous region, e.g. a collection of cells made through extract_region? (2) Could this be implemented such that, if you were using Zeus, the universal fields "x-velocity", "y-velocity" and "z-velocity" were over-written by these new definitions? (If that was done throughout yt, then other references to Zeus in routines like DivV would have to be removed. This is probably still easier than adding in new ones everywhere). (3) Is it possible for me to do that on my own local version? e.g. can I put field definitions in a folder to over-write the default yt ones? Thanks~ Elizabeth

2 2

Proposal: deprecate lazy_reader
by Matthew Turk Feb. 23, 2012

Feb. 23, 2012

Hi all, Long ago, I added the "lazy_reader" option to profiles. This then showed up also in DerivedQuantities. In almost all cases used, this is set to True, and I believe that the only case when it cannot be is when calculating energy. lazy_reader = True means data is read on demand, grid by grid (mediated by preloading for those formats where multiple grids might be in a single file) and False means that all of the data necessary for the reduction is loaded immediately into memory. I'd like to propose we deprecate it as an option, and force all profiles to be "lazy". I would be implementing this in the geometry handling branch. In the geometry handling branch, I have implemented a chunking scheme which is more flexible; in this way subsets of data objects can be loaded. By moving to this scheme, the chunks can either be the *entire* data object (i.e., the non-lazy method from before) or subsets divided in several different ways. For example, you could do grid-by-grid, you could do gridfile-by-gridfile, or whatever type of scheme you like. I believe that getting rid of lazy and moving to just always operating on chunks is the best way to accomplish this. [+-][01] on removing lazy_reader as optional? -Matt

5 6

Parallelism and iteration
by Matthew Turk Feb. 23, 2012

Feb. 23, 2012

Hi all, As some of you may know, I'm working on the geometry subsystem. Right now I'm in the process of constructing incremental reads ("chunking") off disk, with intermediate reduction operations. In the past, this was usually handled by individual tasks, such as profiles, in a pretty gross and manual way. In fixing this up, I noticed we have three different parallel iteration methods. This question is probably best directed at Sam or Britton, but others might also have an opinion. Is there any reason why, from within a subclass of ParallelAnalysisInterface, I wouldn't just use parallel_objects instead of the GridObjectIterator or GridIterator? As a related question, can I get rid of the warning for when running parallel_objects on a single processor? (The geometry handling is coming along nicely, I think. Right now reading from objects works well (and is faster) and the abstraction of grids as opposed to octs or particles has worked well. It's a ways from being feature complete, but if you want to play with it it's all in bb://MatthewTurk/yt-refactor .) -Matt

2 1

call for testing: Orion particle support
by j s oishi Feb. 22, 2012

Feb. 22, 2012

Hi All Orion Users, I've updated my pull request (#99, https://bitbucket.org/yt_analysis/yt/pull-request/99/handle-case-of-orion-d…) to fix Orion particle support for the cases of Star, Sink, and no particle files. Could some of you please test this code and see if it works? Once this is done, someone can accept the pull request. thanks, J.S.

2 4

RFC: introduce {length, mass}_units (similar to time_units)
by Kacper Kowalik Feb. 22, 2012

Feb. 22, 2012

Hi, I've been hit by a part of the code (plot_modifications:get_smallest_appropriate_unit) that assumed my data knows what a Mpc, or even Rsun is. While it's trivial to workaround by using try/except, that made wonder why time_units get a special treatment in StaticOutput? Would it be beneficial to introduce corresponding {lenght,mass}_units? Cheers, Kacper

4 7

Call for homepage viz
by Matthew Turk Feb. 15, 2012

Feb. 15, 2012

Hi all, I'd like to update the viz on the yt-project.org home page. If you have some yt-created visualization you're particularly proud of, please drop me a line off-list with a reasonably sized image and appropriate credit, and if it fits the feel of the site, I'll put it up in the main rotation. A few of the existing images will get rotated out, too. Thanks, Matt

1 0

Parallelization & the ART frontend
by Christopher Moody Feb. 9, 2012

Feb. 9, 2012

Hi guys, I've been working hard on the ART frontend. Lately, I'm to the point where I'm playing around with more complex datasets that are taking much longer to project - so I'd really like to to start using the parallelization engines. I've tried Sam's workshop parallelization demos, and they all work. But launching with the ART frontend ( http://paste.yt-project.org/show/2152/) spawns many independent processes which evidently are not actually splitting the projection job, but still taking up lots of processors. My mpi installation works: yt : [INFO ] 2012-02-07 18:12:28,207 Global parallel computation enabled: 0 / 8yt : [INFO ] 2012-02-07 18:12:28,207 Global parallel computation enabled: 2 / 8 yt : [INFO ] 2012-02-07 18:12:28,208 Global parallel computation enabled: 1 / 8 yt : [INFO ] 2012-02-07 18:12:28,208 Global parallel computation enabled: 6 / 8 yt : [INFO ] 2012-02-07 18:12:28,208 Global parallel computation enabled: 3 / 8 yt : [INFO ] 2012-02-07 18:12:28,208 Global parallel computation enabled: 4 / 8 yt : [INFO ] 2012-02-07 18:12:28,208 Global parallel computation enabled: 5 / 8yt : [INFO ] 2012-02-07 18:12:28,209 Global parallel computation enabled: 7 / 8 But the script is just run 8 times, not any faster. What am I missing here? Many thanks! chris

2 7