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