[Yt-dev] Radius Calculation

[Matthew Turk]

Hi guys,

I did a bunch of testing on this today and satisfied myself it was
correct.  I also fixed a regression I introduced, clarified the
setting of field parameters, and then made AMRSpheres periodic.  This
was committed in r1267.  Code review is welcome!

-Matt

On Fri, Apr 17, 2009 at 9:10 PM, Matthew Turk <matthewturk at gmail.com> wrote:
> Okay, so suddenly it occurs to me -- if your domain width is 1.0, then
> this is in fact just a subset of your domain, and it could be correct.
>  I forgot about parallel root grid IO.  What happens if you plot a
> slice through the center of the domain of the new field?
>
> pc = PlotCollection(pf, center=[0.5,0.5,0.5])
> pc.add_slice("NewRadius", 0)
> pc.add_slice("Radius", 0)
> pc.save("new")
>
> That should show us if it's actually giving the right results or not.
> Apologies; I completely forgot about decomposed root grids when I
> wrote that script.
>
> On Fri, Apr 17, 2009 at 9:01 PM, Matthew Turk <matthewturk at gmail.com> wrote:
>> Hm, this looks wrong to me.  Clearly something I have not considered
>> -- over here myi mages look like the attached.  What's the value for
>> the domain width?  This is confusing to me, as I really thought it
>> would work!
>>
>> Can you try with a projection?  It should work, right through the
>> domain, and it should tell us what we're looking for.  Or, jsut
>> straight-up slices of NewRadius.  I'll try some stuff too.
>>
>> On Fri, Apr 17, 2009 at 8:48 PM, Britton Smith <brittonsmith at gmail.com> wrote:
>>> I just tried your script with some data I have here.  I've attached the
>>> images.  I can't really tell if this is right or not.  I think another way
>>> to test this would be to do a projection with a periodic sphere as the
>>> source object.  If you throw a sphere that overlaps the domain boundary to a
>>> projection, it should be quite clear if it's doing what it's supposed to.
>>> Does that sound right?
>>>
>>> Britton
>>>
>>> On Fri, Apr 17, 2009 at 8:17 PM, Matthew Turk <matthewturk at gmail.com> wrote:
>>>>
>>>> Oops, just as a note, there's an error in the paste -- the second plot
>>>> of "NewRadius" should be of "Radius".
>>>>
>>>> Thanks guys...
>>>>
>>>> -Matt
>>>>
>>>> On Fri, Apr 17, 2009 at 7:16 PM, Matthew Turk <matthewturk at gmail.com>
>>>> wrote:
>>>> > 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 at 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-new-universal-function
>>>> >>
>>>> >> 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
>>>> >>
>>>> >
>>>> _______________________________________________
>>>> 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
>>>
>>>
>>
>