On 7 Oct 2013, at 21:55, Ian Hinder <[email protected]> wrote:

> 
> On 7 Oct 2013, at 19:28, Roland Haas <[email protected]> wrote:
> 
>> Ps4 mode decomposition:
>> * Jim Healy sees some odd behaviour in 2,1 mode of psi4 when using
>> Multipole thorn but not when using old GT mode decomposition thorn
>> * mode that is expected to be zero is not and leads to BH kick
>> * run uses rotatingsymmetry180 and reflectionsymmetry (in z direction)
>> * no physics change in multipole since Oersted release
>> * suggest running in full mode
>> * will work with Jim for a while then possibly take issue to mailing
>> list again
> 
> The other aspect was that this worked OK with the Oersted release, but fails 
> with the Gauss release.  Since there were no significant changes in Multipole 
> or AEILocalInterp between those two releases, the cause probably lies 
> elsewhere (e.g. Carpet).
> 
> Jim, could you try running without symmetries, to see if the problem is 
> related to symmetries?  Also, you could try setting the parameter 
> CarpetLib::interpolate_from_buffer_zones = yes in the Gauss version?  The 
> default changed from "yes" to "no" between these two releases. 


According to Jim, the problem was not in fact newly-introduced; the same 
results are obtained with earlier releases. There must have been some confusion 
in the telecon.  I have investigated, and the culprit appears to be in one of 
the integration methods used by Multipole.  Multipole interpolates data onto a 
sphere at regular intervals in theta and phi, and then performs a 2D numerical 
quadrature on these points.  It has several integration methods available, 
which give results of various orders of accuracy.  The default (first order 
accurate) integration method, labelled "midpoint" [** see below], is 
implemented in an unusual way.  In the following, I will discuss 1D 
integration, but the generalisation to 2D is straightforward.

The natural first order accurate approximation of an integral is given by 
summing the integrand at all the points from a to b-dx and multiplying by dx, 
i.e. by adding up the areas of the rectangles.  This is the analogue of the 
Riemann sum integral definition.  However, in Multipole, the sum runs from a to 
b; i.e. it "over counts" by one.  The effect of this is to introduce an 
additional first order error.  So the resulting error is still first order 
(which is what was originally tested when Multipole was developed), but it is 
not the most natural definition.  Additionally, this method does not have the 
property that a function antisymmetric about (a+b)/2 integrates numerically 
exactly to zero, though the integral will converge to zero as O(h)), just as 
the integral of any other function converges to the exact result as O(h).

In Jim's case, increasing the number of angular points in the phi direction 
should make the odd-m modes converge to zero if the underlying data is 
Pi-symmetric.  Alternatively, the other integration methods could be used, as 
they satisfy the symmetry property.

Typically, we set the number of points on the sphere to be so large that these 
issues should be negligible in comparison to other (accumulating) sources of 
error, but Jim's experience reminds us that this should be checked.  I have 
tended to worry that higher order integration methods might lead to more noise, 
but I haven't checked this.  It might be wise to use the higher order methods.

** Note on naming: the Wikipedia page on the "rectangle" or "midpoint" method, 
http://en.wikipedia.org/wiki/Rectangle_method, is very confused I think.  The 
expression given there is not a sum over the midpoints of the intervals, but is 
a sum over the left points of each interval.  The "Error" section then goes on 
to claim that the error is O(h^2).  This is true if you sum over the midpoints, 
but NOT if you sum over the left points (you can check this using a Taylor 
expansion).  If you sum over the left points, you get an O(h) accurate 
approximation of the integral.  I suspect that the first order accurate method 
in Multipole was named "midpoint" after a lazy wikipedia reference (probably by 
me).

Should we do anything about this?  The default "midpoint" method in Multipole 
is not only very badly named, it also does not have the desired behaviour of 
integrating an antisymmetric integrand exactly to zero.  

I think we should deprecate this integration method, with a warning that the 
default is going to change to "simpson" in future.  We could optionally 
implement a "rectangle" (is this the right name?) method which is the same as 
the current first order method but excluding the endpoints.

-- 
Ian Hinder
http://numrel.aei.mpg.de/people/hinder

Attachment: signature.asc
Description: Message signed with OpenPGP using GPGMail

_______________________________________________
Users mailing list
[email protected]
http://lists.einsteintoolkit.org/mailman/listinfo/users

Reply via email to