On 21 Oct 2010, at 20:58, Bernard Kelly wrote:
Hi all.
I'm trying some GW extraction from a Carpet/McLachlan black-hole binary evolution. I thought the combination of thorns WeylScal4 and Multipole would give me something reasonable, but I'm having trouble with it.
I'm attaching a parameter file for my run. As you can see, the extraction is supposed to be over a set of coordinate spheres at radii 45, 50 55, 60 65, & 70.
What I get is NaN in the WFs at all timesteps after t=0. I'm also outputting 2D data (though not as frequently), and Psi4 seems well- behaved. So what's going wrong?
I have noted that though the routine "Multipole_Calc" should happen after "calc_np" according to its schedule.ccl, it appears -before- GROUP WeylScal4_Calculate in the run's initial printing of scheduled routines.
Could you post the schedule output? According to the copy of the code that I have from the release version, Multipole_Calc is called in ANALYSIS and WeylScal4_Calculate is called in MoL_PseudoEvolution, which is called in EVOL (and some other places). The calc_np is a relic from when both were called in ANALYSIS.
Suggestions appreciated,
I tried to run your parameter file to test, but unfortunately I have some work-in-progress hacks to my local copy of WeylScal4, so it didn't run. Nevertheless, I think I know what the problem is. The spatial triad used by WeylScal4 is a coordinate triad (d/dph, d/dr, d/ dth), suitably orthonormalised. This triad is singular on the z axis, and the way it is computed there leads to a NaN. When you compute spherical harmonic modes, points on the z axis (th = 0) are multiplied by sin(th) = 0, so the value on the axis is never used, but numerically the NaN still causes problems.
There are several approaches to this:
1. Replace the triad at that point with a different one. This way, Psi4 on th = 0 is still "Psi4", just with respect to a different tetrad at that point. This is probably the best approach, and I have it implemented in my local copy. I can commit this to the development version.
2. Use an open integration formula on (0,Pi) on the 2-sphere when doing the mode decomposition. This means that the point at th = 0 is not used. Multipole does not do this, but it wouldn't be too hard to modify it.
3. Do a numerical trick to make sure that the point at th = 0 doesn't have a NaN. The particular trick I usually use is to offset the coordinates used to construct the spatial triad by 10^-15 in one direction.
WeylScal4::offset = 1e-15
This is what I did before I implemented the triad change on the axis.
What do other people do to get around this problem? Would some more physical tetrad (e.g. quasi-Kinnersley?) help?
P.S. Is this appropriate to einsteintoolkit, or Cactus at large? I'm using the ET_2010_06 release ...
WeylScal4 and Multipole are both part of the Einstein Toolkit (in the EinsteinAnalysis arrangement), so discussion is probably best on the ET list. I have CC'd this message to that list.
Other comments about your parameter file:
* You have asked for the mode decomposition of Psi4r and Psi4i individually, and you have not set the spin weight parameter, so it will default to 0, which is not what you want. Try
Multipole::variables = "WeylScal4::Psi4r{sw=-2 cmplx='WeylScal4::Psi4i' name='Psi4'}"
instead. This tells Multipole to use spin weight -2 for Psi4, and to treat Psi4r and Psi4i as the real and imaginary parts of a complex number. This means you don't have to recombine the modes of the real and imaginary parts by hand. WeylScal4 and Multipole could be modified to work with Cactus complex grid functions, but I don't have any experience with using them, and Kranc does not currently support them.
* Your time_refinement_factors parameter is quite aggressive - you are using time refinement only for finest four levels. You can probably get away with time refinement for the last 6, or more depending on your finite differencing order and value of eta.
* I recommend using the parameters
TimerReport::out_every = 32 TimerReport::n_top_timers = 40
which will give you a quick overview of the parts of the code which are taking the most time.
Apart from those points, the parameter file looks fine.