Hi guys, Can somebody review this? I chatted with Jeff, and I believe it is correct for periodic boxes, but I need somebody to look at it before I commit it, as it is a rather crucial aspect of yt. This changes the way radius is calculated and should now work with periodic boundary conditions. Here's the new field definition: def _NewRadius(field, data): center = data.get_field_parameter("center") DW = data.pf["DomainRightEdge"] - data.pf["DomainLeftEdge"] radius = na.zeros(data["x"].shape, dtype='float64') for i, ax in enumerate('xyz'): r = na.abs(data[ax] - center[i]) radius += na.minimum(r, na.abs(DW[i]-r))**2.0 na.sqrt(radius, radius) return radius the final argument in na.sqrt is the 'out' argument, to avoid temporary arrays. I've committed a sample script here: http://paste.enzotools.org/show/101/ which loads up a data dump, calculates the radius from [0.1, 0.1, 0.1] at all points in the root grid, then plots the output. It gives correct min & max values for the radius. If I can get somebody else to sign off on this, I will commit the necessary changes to _Radius and _ParticleRadius, as well as grid & point selection for spheres. We'll then have a fully-periodic profiler, which will close #165. -Matt On Thu, Apr 16, 2009 at 11:39 AM, Matthew Turk <matthewturk@gmail.com> wrote:
Hi guys,
I've been exploring the idea of changing radius to be correct for periodic boxes. Right now it is incorrect; each component (x,y,z) should not have a distance greater than 0.5 * domain_size. The easiest way to do this would be:
rx = min( abs( x - center_x ), abs( x - center_x - domain_x) )
I think the best way to do this is to set up a NumPy ufunc
http://docs.scipy.org/doc/numpy/user/c-info.beyond-basics.html#creating-a-ne...
that accepts three arrays, along with the center and the domain width, and then returns the radius. What do you all think? Alternatively, I'm thinking maybe just a regular function that gets the arrays and returns one would be better; the ufunc machinery is a bit complicated and I might get confused. Once I come up with it, will somebody be able to look over my work?
-Matt