[FLASH-USERS] Indexing Bug in Multipole_experimental

[Norbert Flocke]

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
>