[yt-dev] Matching points to grids

[Matthew Turk]

Hi John,

On Sun, Nov 18, 2012 at 6:28 PM, John ZuHone <jzuhone at gmail.com> wrote:
> 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 completely agree, this would be very, very useful, particularly for
initial condition generation on the whole.  And yes, you're right,
we'd need to generalize a fast version of find_point.

Looking at the code, there are a few optimizations that could be put
in place -- for instance, you could get rid of allocations by doing
in-place logical_and operations.  Like:

cond = np.ones(pf.h.num_grids, dtype="bool")
for ax, p in enumerate( [x,y,z] ):
    np.logical_and( pf.h.grid_left_edge[:,ax] <= p, cond, cond)
    np.logical_and( pf.h.grid_right_edge[:,ax] >= p, cond, cond)

This still allocates a temporary boolean array for the conditional
inside logical_and, and there's probably a way around that, but I
can't recall it.

I would like most to see -- in the long run -- a Cython data structure
that represents the grid/oct structure that can simply be walked.
This would speed up nearly all of our geometric selection.  The best
place for it is 3.0 I believe.

>
> 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.

I love this idea.

>
> 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.

This is also quite nice.  I've written something like this for an Enzo
output (to add tracer particles) but it was not general.

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

I think the path you are on is a good one at the moment.  The only
thing I can see that we'd want to swap out would be the actual
gridding, not the particle generation or the writing/loading
overrides.

As another note, the idea of a "sidecar" file that can be loaded up
along side a dataset is *really* compelling.  Imagine writing out
fields separately from the simulation code and using them to override
the sim code itself.  This could be used for things like tagging
individual cells, for instance.  And the machinery you're developing
would address that nicely.

-Matt

PS Again sorry for all the additional emails that got sent out ...
embarrassing, although I still blame the GMail Offline extension!

>
> 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
>
>
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
>