[FLASH-USERS] Indexing Bug in Multipole_experimental

[James Guillochon]

Hi Norbert,

Even with the latest suggested change, I have found that the value of the
potential calculated is quite inaccurate if the grid cell center happens to
lie less than grid_cell_size/inner_zone_size distance away from the center
of mass. I think there still must be a bug in how the potential is
calculated when r --> 0. I've been pouring over the code for a few days now,
and I haven't been able to determine exactly what to change to improve the
situation...would you be so kind to look at the code and see if there's
anything that can be done to improve the potential calculation right next to
the center of mass?

(FYI, I have tried increasing the inner_zone_size, but in my calculation the
center of mass moves quite slowly across the grid, and inevitably there is
time step in the calculation where r << grid_cell_size/inner_zone_size).

Thanks!
- James

On Wed, Sep 21, 2011 at 8:33 AM, Norbert Flocke
<flocke at flash.uchicago.edu>wrote:

> Hi James and all FLASH users,
>
> Yesterdays comment of James Guillochon about an indexing problem arising
> for small values of the variable 'rscaled' in the gr_mpolePot* and
> gr_mpoleMom* routines, even when the rscaled values are not near the
> computational zero, had a red flag rising in my head. After having
> a closer look at the code, the fix is not the one suggested by James
> (sorry), but the following needs to be done:
>
>  In the routine 'gr_mpoleSetInnerZoneGrid' go to the very end and
>  replace the section ...Create the inner zone radii locator... by the
>  following:
>
> !
> !
> !     ...Create the inner zone radii locator arrays to help search through
> !        the inner zone grid. This is done on all processors.
> !
> !
>  allocate (inner_zone_Qlower (0:inner_zone_size))   ! lower index 0
>
>  inner_zone_Qlower = 0
>
>  Qstart = 1
>  do n = 0,inner_zone_size                           ! start from 0
>     do Q = Qstart,inner_zone_qmax
>        dr_unit = int (inner_zone_radii (Q))
>        if (dr_unit > n) then
>            exit
>        else if (dr_unit == n) then
>            inner_zone_Qlower (n) = Q
>            Qstart = Q + 1
>            exit
>        end if
>     end do
>  end do
>
>  allocate (inner_zone_Qupper (0:inner_zone_size))   ! lower index 0
>
>  inner_zone_Qupper = 0
>
>  Qstart = inner_zone_qmax
>  do n = inner_zone_size,0,-1                        ! end at 0
>     do Q = Qstart,1,-1
>        dr_unit = int (inner_zone_radii (Q))
>        if (dr_unit < n) then
>            exit
>        else if (dr_unit == n) then
>            inner_zone_Qupper (n) = Q
>            Qstart = Q - 1
>            exit
>        end if
>     end do
>  end do
>
> The changes have been marked with a ! comment. As you can see the
> 'inner_zone_Qlower' and 'inner_Zone_Qupper' arrays are reindexed
> and the search is extended to include the '0' index. The routines
> gr_mpolePot* and gr_mpoleMom* need not to be changed.
>
> Thank you very much.
> Norbert
>
>
>
> On Tue, 20 Sep 2011, James Guillochon wrote:
>
>  Hi Norbert, just to clarify, rscaled is not zero in this case, it's just a
>> small number. The issue is that the ceiling function only applies to the
>> first product in the expression rscaled*inner_zone_grid_inv, and then the
>> int() rounds the second product to the nearest whole number.
>>
>> On Tue, Sep 20, 2011 at 10:09 AM, Norbert Flocke
>> <flocke at flash.uchicago.edu>**wrote:
>>
>>  Dear James Guillochon,
>>>
>>> You are right. This is a computer representability issue for very small
>>> numbers. In case r falls below the smallest representable number, then
>>> rscaled is zero and the ceiling functions returns 0 as well. Thank you
>>> very much for reporting this issue. It will be fixed right away.
>>>
>>> Regards,
>>> Norbert Flocke
>>>
>>>
>>>
>>> On Tue, 20 Sep 2011, James Guillochon wrote:
>>>
>>>  Hi, I've discovered a small bug in Multipole_experimental. In a couple
>>>
>>>> places in gr_mpolePot* and gr_mpoleMom* there is the following set of
>>>> lines:
>>>>
>>>> rscaled = r * dr_inner_zone_inv
>>>>
>>>>  dr_unit = int (ceiling (rscaled * inner_zone_grid_inv) *
>>>>> inner_zone_grid)
>>>>> Q_lower = inner_zone_Qlower (dr_unit)
>>>>>
>>>>>
>>>>>  Occasionally, if r is small enough, dr_unit will be 0, which causes an
>>>> out
>>>> of bounds error when referencing inner_zone_Qlower. The fix is to
>>>> replace
>>>> the dr_unit assignments with the following:
>>>>
>>>> dr_unit = max(int (ceiling (rscaled * inner_zone_grid_inv) *
>>>>
>>>>  inner_zone_grid), 1)
>>>>>
>>>>>
>>>>>  Cheers,
>>>> - James
>>>>
>>>> --
>>>> James Guillochon
>>>> Department of Astronomy & Astrophysics
>>>> University of California, Santa Cruz
>>>> jfg at ucolick.org
>>>>
>>>>
>>>>
>>
>> --
>> James Guillochon
>> Department of Astronomy & Astrophysics
>> University of California, Santa Cruz
>> jfg at ucolick.org
>>
>>


-- 
James Guillochon
Department of Astronomy & Astrophysics
University of California, Santa Cruz
jfg at ucolick.org
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20110930/1c2613bd/attachment.htm>