On 9 Oct 2013, at 22:48, Erik Schnetter <[email protected]> wrote:
> On 2013-10-09, at 16:28 , Ian Hinder <[email protected]> wrote: > >> On 9 Oct 2013, at 20:14, Erik Schnetter <[email protected]> wrote: >> >>> On 2013-10-09, at 13:16 , Ian Hinder <[email protected]> wrote: >>> >>>> 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. >>> >>> The domain in the phi direction is periodic. Introducing an arbitrary >>> boundary in the phi direction does not make sense. The only sensible >>> integration method is to take the sum over all points in the phi direction >>> (and to divide by N). >> >> In our case, we interpolate from 0 to 2 Pi, so we need to exclude the last >> (repeated) point. Come to think of it, this is fairly dumb for the phi >> direction, but it makes the code simpler to deal with theta and phi >> uniformly. Maybe this should be changed (at the cost of making all >> waveforms computed previously unreproducible). > > The interval is [0, 2pi); see line 295 in setup.cc. Where is this file? It's not in Multipole. Are you thinking of SphericalSurface? Multipole does not use that thorn; it aims to be independent. > If you have N interior points (excluding ghosts), then the first point is a > phi=0, and point N (which is the first ghost point past the domain) is at > phi=2pi. The notation in Multipole is that N is the number of intervals, not the number of points. There are N+1 points, and in the phi direction, the last one and the first one are the same. > If you integrate while assuming a boundary, then the boundary points are 0 > and N. For example, applying the midpoint rule leads to > > s := (1/N) sum[i=0..N-1] (1/2) (x[i] + x[i+1]) > > which is equivalent to > > s := (1/N) sum[i=0..N-1] x[i] > > after replacing x[N] with x[0]. Of course, if the domain is not periodic, > then this doesn't work. > >>> If you do not want to take this as granted, then consider the following: >>> - choose are arbitrary (!) integration method in the phi direction >>> - perform this integration N times, moving the integration boundary by 1 >>> point for each step. this is possible since the domain is symmetric. each >>> point then becomes the boundary for one of these steps. >>> - take the average of these N integrations >>> - by linearity, these N integrations can be collapsed to a single >>> integration >>> - by symmetry, the coefficient for each grid point for this collapsed >>> integration has to be the same >>> - qed >> >> This just shows that there exists an integration method which has all the >> weights the same, which is already known. I don't think it says anything >> about the individual methods. > > Since we started from an arbitrary method, and since the choice of boundary > does not influence convergence in a periodic domain, this shows that > equal-weights is at least as good as any method. Hence it is no worst than > the best method. > >>> Put differently: In a periodic grid, no point is special. Hence, the >>> integration weight for each point must be the same. >>> >>> This corresponds to a spectral method (using a Fourier basis), with >>> exponential convergence. >>> >>> We should remove all other phi integration methods. >> >> It looks like you know much more about this than I do. I didn't know about >> this result for quadrature of periodic functions, but now that I think about >> it some more, it makes sense. The integration methods in Multipole are the >> classical quadrature methods for \int_a^b f(x) dx with no additional >> information about f. >> >> My worry about using a highly accurate (i.e. high order) method with a small >> number of points is that noise at some of these small number of points could >> significantly affect the result. Do you know anything about this? Am I >> worrying for nothing? > > > This is correct, the noise will influence the results. However, if you are > worried about noise, then you should not be looking for high-order methods; > instead, a low-order method with many points will lead to smoother results. > Thus, you would still choose a large number of points in the phi direction to > avoid the noise, ignoring any arguments about spectral convergence. > > In the theta direction, the noise can be problematic, since the integration > weights become very oscillatory for large N. It may be better there to > restrict ourselves to (say) 8th order. > > -erik > > -- > Erik Schnetter <[email protected]> > http://www.perimeterinstitute.ca/personal/eschnetter/ > > My email is as private as my paper mail. I therefore support encrypting > and signing email messages. Get my PGP key from http://pgp.mit.edu/. > -- Ian Hinder http://numrel.aei.mpg.de/people/hinder
signature.asc
Description: Message signed with OpenPGP using GPGMail
_______________________________________________ Users mailing list [email protected] http://lists.einsteintoolkit.org/mailman/listinfo/users
