[yt-dev] Matching points to grids

Sun Nov 18 15:28:29 PST 2012

Hi Matt,

I suppose I should get into why I've been working on this. I had a need to initialize some tracer particles in a FLASH simulation that's already been running. I decided to use yt to open up the checkpoint file and generate the tracer particles based on the gas density, and then I was going to write them back to the checkpoint file in the format FLASH expects.

It then occurred to me that such particle initialization could be useful for yt in general, particularly for initial conditions generation as you and others have been working on. So I have been working on some classes for generating particle fields on the fly, with positions based on some criterion such as density, a lattice, or any general criterion. For some distributions, the grids the particles belong to would not be known a priori, so it seemed necessary to write a faster implementation of find_point and generalize it to a number of points.

I was thinking about two applications in particular:

1) One could generate a grid structure with some fields, generate particle positions, generate other particle fields, and then write everything to disk using GDF.

2) An existing pf has no particles, but the user has some other criterion that would be used to generate a set of particle fields that they are interested in. These fields could then be "loaded in" as any others would.

Any thoughts on this? I've already written most of the machinery for this in the last few days.

John 

John ZuHone
Laboratory for High-Energy Astrophysics
NASA/Goddard Space Flight Center
8800 Greenbelt Rd., Code 662
Greenbelt, MD 20771
(w) 301-286-2531
(m) 773-758-0172
jzuhone at gmail.com
john.zuhone at nasa.gov

On Nov 18, 2012, at 5:56 PM, Matthew Turk <matthewturk at gmail.com> wrote:

> Hi John,
> 
> I think this looks pretty good for now, certainly for 2.x.
> 
> In 3.0 I want to eventually transition FLASH to the Octree geometry structure, which would mean that particle location could all occur directly inside Cython.  We'd also have a number of much faster particle selection routines, which could eventually make it not necessary to directly sort or select particles too strongly in advance, as long as a coarse index could be built.
> 
> I say push it and we can move on from there.  If you want, one thing you could absolutely do right now is construct a Cython routine that accepted:
> 
> x,y,z
> left_edges, right_edges, parent_index
> 
> and then construct an octree from that, assigning particles as they came.
> 
> -Matt
> 
> On Saturday, November 17, 2012, John ZuHone wrote:
>> By the way, I certainly think that there is much room for improvement, so anyone that would like to take a look at it and make suggestions is most welcome. 
>> 
>> On Nov 17, 2012, at 6:09 PM, Christopher Moody <cemoody at ucsc.edu> wrote:
>> 
>>> Hi John,
>>> That speed is similar to what I found. Brute force collision check would make (3x2x7000x10^6) ~ 10^10 floating point operations, which should take a few seconds on a modern single core (3Ghz*10seconds). Making an actually intelligent depth-first search should drop the number of grids to check for every particle by an order of magnitude. So we should really be able to beat a 4-minute gridding time, we just need to figure out what the bottleneck is. And of course this is applicable to particles in yt-3.0.
>>> 
>>> chris
>>> 
>>> 
>>> On Sat, Nov 17, 2012 at 6:30 AM, John ZuHone <jzuhone at gmail.com> wrote:
>>> Chris,
>>> 
>>> Since last night I've developed a working depth-first approach which is *reasonably* fast, but I think it depends on your definition of reasonable. :) It just did a FLASH dataset with ~7000 blocks and 10^6 points in about 4 minutes. 
>>> 
>>> I'm going to try it out on a few more datasets, and I'm hoping to issue a PR for it today or tomorrow. 
>>> 
>>> Best,
>>> 
>>> John
>>> 
>>> On Nov 16, 2012, at 6:21 PM, Christopher Moody <cemoody at ucsc.edu> wrote:
>>> 
>>>> Hi John, 
>>>> I tried writing a cython routine that matched particles to grids a while back. It works in the opposite order as a recursive octree-like gridding. It starts with the highest-level most-refined grids, recording the grid ID that each particle belongs to. It skips previously gridded particles as it progresses to coarser grids. It's likely not faster than depth-first approach.
>>>> 
>>>> It could be helpful, but it's terribly slow, despite several attempts to speed it up. 
>>>> Check out  assign_particles_to_cell_lists() in yt-hg/yt/utilities/lib/CICDeposit.pyx
>>>> 
>>>> chris  
>>>> 
>>>> 
>>>> On Fri, Nov 16, 2012 at 11:30 AM, Matthew Turk <matthewturk at gmail.com> wrote:
>>>> Hi John,
>>>> 
>>>> I think that would work fine for Enzo, we'd just need to iterate over
>>>> root grids.
>>>> 
>>>> -Matt
>>>> 
>>>> On Fri, Nov 16, 2012 at 2:28 PM, John ZuHone <jzuhone at gmail.com> wrote:
>>>> > Hi Matt,
>>>> >
>>>> > Would it work just to start at the top-level grids and then work down the tree, checking the children as you go? I know this would work for FLASH, but I'm not certain it would work for Enzo.
>>>> >
>>>> > John
>>>> >
>>>> > On Nov 16, 2012, at 2:14 PM, Matthew Turk <matthewturk at gmail.com> wrote:
>>>> >
>>>> >> Hi John,
>>>> >>
>>>> >> That's the ideal way, but we don't yet have the code to do that.  I
>>>> >> would like to do that.  The fastest method would be to create a new
>>>> >> version of find_values_at_points (which accepts an array) and make itj
>>>> >> ust return grid indices.  This would do the NxM collision check of
>>>> >> particles in grids.
>>>> >>
>>>> >> Writing an octree-aware particle finder would be a good idea.
>>>> >> Necessary, even...
>>>> >>
>>>> >> -Matt
>>>> >>
>>>> >> On Fri, Nov 16, 2012 at 1:07 PM, John ZuHone <jzuhone at gmail.com> wrote:
>>>> >>> I did know about that, but given, say, a million points, would that be
>>>> >>> *fast*? I need to look at what it does I guess.
>>>> >>>
>>>> >>> In FLASH, what'd we'd normally do is start at the top level and traverse
>>>> >>> down the octree until we found the block.
>>>> >>>
>>>> >>> On Nov 16, 2012, at 1:02 PM, Nathan Goldbaum <nathan12343 at gmail.com> wrote:
>>>> >>>
>>>> >>> Hi John,
>>>> >>>
>>>> >>> Yes, you can use pf.h.find_point() to do that.  There's some documentation
>>>> >>> on it here:
>>>> >>> http://yt-project.org/doc/analyzing/low_level_inspection.html#finding-data-at-fixed-points
>>>> >>>
>>>> >>> -Nathan
>>>> >>>
>>>> >>> On Nov 16, 2012, at 10:00 AM, John ZuHone wrote:
>>>> >>>
>>>> >>> Hi all,
>>>> >>>
>>>> >>> Do we have anything in yt right now that, given a set of points, *quickly*
>>>> >>> finds out which grids they belong to?
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-dev-spacepope.org/attachments/20121118/434a4be6/attachment.htm>