[FLASH-USERS] Indexing Bug in Multipole_experimental
Norbert Flocke
flocke at flash.uchicago.edu
Thu Oct 6 10:44:30 EDT 2011
All FLASH users,
The apparent indexing bug in the Multipole_experimental solver has
been resolved. There is no need to fix anything in the code. The
bug was caused by addition of an extra routine to the solver on a
local copy of FLASH. For anyone interested in the details, please
consult the relevant emails exchanged (see below, backwards
chronological order).
Thank you very much,
Norbert Flocke
-----------------------------------------------------------------------------------
On Mon, Oct 3, 2011 at 9:00 AM, Norbert Flocke
<flocke at flash.uchicago.edu>wrote:
James, I now know what is happening. I will try to explain.
The reason you are getting errors is the following:
1) When the inner radial zone is defined and set up by the
Multipole solver, it constructs additional arrays to help
search in a fast way to which Q (this is a radial bin in
the inner zone) the current {x,y,z} coordinate belongs to.
However, and this is the important point, these additional
arrays (inner_zone_Qupper + inner_zone_Qlower) are constructed
such that they assume that the only possible {x,y,z} points
correspond to the mid points of cells. The code does NOT
know that you are about to calculate the potential at arbitrary
points. With other words, the code assumes that all possible
{x,y,z} fall into one of the inner zone radial bins. If you
request the code to calculate into which bin your arbitrary
{x,y,z} triple will fall, it must fail, unless by coincidence it
actually lies in one of the bins. This is no problem for the
outer zone bins, as their boundaries are touching each other
in space and henceforth your {x,y,z} must per force be in
one of the outer zone bins.
2) In order to get the code do what you want in the inner zone,
you have to look at the formula that is used to calculate the
potential for the 'normal' case when the {x,y,z} is in a bin.
It reads (the || bars denote one bin):
|Q-1| | Q | | Q+1 |
^
{x,y,z}
Potential{x,y,z} = R dot MI(Q+1) + I dot [MR(Q) - m*R]
In this formula, R and I are the {x,y,z} R and I moments,
MI and MR are is the cummulative sum of I and R moments and
m*R is the correction you have to substract because your
{x,y,z} is in the Q-th bin. So far so good. Now, if your {x,y,z}
is outside of any of these bins, then we have:
|Q-1| | Q | | Q+1 |
^
{x,y,z}
Potential{x,y,z} = R dot MI(Q+1) + I dot MR(Q)
Notice, no correction is necessary here, but you must also
find the right bins Q and Q+1 (which the code is NOT equipped
to do). A special situation arises, if your {x,y,z} is between
the first inner zone bin Q=1 and the center of mass:
CM | Q=1 |
^
{x,y,z}
Then the potential is:
Potential{x,y,z} = R dot MI(Q=1)
I hope this clarifies the situation. Unfortunately it means more coding
to catch the situations described under 2).
Cheers,
Norbert
---------------------------------------------------------------------------
On Fri, 30 Sep 2011, James Guillochon wrote:
Norbert,
I created a custom function, gr_mpoleGradPot, which accepts a coordinate
(x,y, z), and then from that it does a finite difference about that point
using the same method as what's used in gr_mpolePot. The finite difference
length is small enough that I can ignore the self-potential as the density
does not change rapidly. It works everywhere in the grid except right next
to the center of mass. I've attached the function so you can see how it works,
and maybe spot any obvious errors.
I appreciate all the help, even though now it seems apparent that the error
is related to using the solver in a non-standard way...
----------------------------------------------------------------------------
On Fri, Sep 30, 2011 at 2:22 PM, Norbert Flocke wrote:
James,
Forgive me, but I still don't understand quite exactly how you calculate
the potential gradient at your point of interest 'o'. The only potentials
that you get from the Multipole solver are potentials associated with
the cell centers (marked 'x' in your picture). To obtain these your current
settings are ok. The Multipole solver cannot give you a potential at the
point 'o'. If you want that then you have to use higher refinement such that
the cell midpoints get closer to 'o'.
-----------------------------------------------------------------------------
On Fri, 30 Sep 2011, James Guillochon wrote:
Hi Norbert,
So I think I may understand the source of the issue. I calculate the
potential gradient using a finite difference template:
phi(x - dx) - phi(x + dx) / (2*dx), where x is *not* necessary a grid
point location, directly from the moments themselves. At the first
time-step, x = 0, and shortly after x = some small number. Thus even
though the inner zone radial gridding is being done based on the
location of the center of mass relative to the grid cell midpoints, I
am trying to calculate the potential gradient at an arbitrary
position that may be closer than any grid cell midpoint to the center
of mass.
Here is a crude drawing showing the situation after step 1:
+-------+-------+-------+-------+
| | | | | * = center of mass
| x | x | x | x | x = cell centers
| | o| | | o = position I want potential gradient
+-------+-------*-------+-------+ -+| = cell boundaries
| | | | |
| x | x | x | x |
| | | | |
+-------+-------+-------+-------+
So in essence, it seems that I have to increase inner_zone_size to
accurately calculate the potential so close to the center of mass, as
based on the grid cell positions alone the inner bin is unlikely to
be small enough to get the potential at the position I want. I've
found that increasing inner_zone_size much above 16 slows things down
tremendously, so alternatively I could just set the innermost bin to
be smaller than the distance between the (*) and the (o) point?
In other words, it appears that the problem is more related to my
meddling than how the code is currently written. Do you think there
is anything I should worry about by setting the innermost radial bin
to be very small?
-----------------------------------------------------------------------------
On Fri, Sep 30, 2011 at 10:44 AM, Norbert Flocke
<flocke at flash.uchicago.edu>******wrote:
Hi James,
I need more information from your run(s). What geometry are you using?
What do you mean exactly by your phrase:
if the grid cell center happens to lie less than
grid_cell_size/inner_zone_size distance away from the
center of mass.
Is 'grid_cell_size/inner_zone_size' a ratio?. The more info you can give
me about your run(s), the easier it is for me to see what is happening.
Regards
Norbert
-----------------------------------------------------------------------------
On Fri, 30 Sep 2011, James Guillochon wrote:
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
-----------------------------------------------------------------------------
More information about the flash-users
mailing list