[yt-users] Definition of the Radius field

[Nathan Goldbaum]

Hi all,

I encountered a bit of odd code deep in the field definition code and I 
need some help to understand why it is the way it is.

The field I'm curious about is the Radius field (defined on line 796 of 
universal_fields.py):

def _Radius(field, data):
     center = data.get_field_parameter("center")
     DW = data.pf.domain_right_edge - data.pf.domain_left_edge
     radius = np.zeros(data["x"].shape, dtype='float64')
     for i, ax in enumerate('xyz'):
         r = np.abs(data[ax] - center[i])
         radius += np.minimum(r, np.abs(DW[i]-r))**2.0
     np.sqrt(radius, radius)
     return radius

In particular, I don't understand this line:

radius += np.minimum(r, np.abs(DW[i]-r))**2.0

Is this supposed to correct for some effect in periodic simulations?

The reason I bring this up is because it caused some very weird behavior 
in a 1D FLASH simulation.  In this simulation, the domain went from x = 
-.01 pc to x = +500 pc.  Since the middle of the domain wasn't lined up 
with the 'center' of the simulation at x = 0, the offending line in the 
radius calculation led to an incorrect calculation of the radius beyond 
x = 250 pc.  In that case, I was able to fix it by changing the 
offending line to

radius += r**2.0

Before I issue a patch, I want to see why it's defined the way it is 
currently and see if we can come up with a workaround.

Cheers,

Nathan