Hi Ian.
Thanks for looking at this so carefully.
On 10/22/10 3:55 AM, Ian Hinder wrote:
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.
Aha. Yes, I changed that back (i.e., put WeylScal4_Calculate back into ANALYSIS), because I wanted both in the same schedule bin. Why did that change? Are variables in MoL_PseudoEvolution available for ANALYSIS automatically (i.e. will they retain their stored values)?
I think part of my problem is that I haven't caught up with the modern set of Cactus bins -- it used to be much simpler. I'm attaching my local copy of WeylScal4's schedule.ccl, as well as the overall scheduling header for the run. The ANALYSIS-bin lines I didn't understand are 415-417 & 620-622.
I've restored the original schedule.ccl now, have recompiled, and am now rerunning. I'll let you know if this helps.
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.
I had a different Carpet grid structure before, and then the same extraction didn't give NaNs. I wonder why?
There are several approaches to this:
- 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.
- 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.
- 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?
I don't think quasi-Kinnersley will be a panacea; theta = 0,pi is still a problem unless you take special precautions. Personally, I think the most "correct" thing to do is your option 2 -- use an open integration formula that avoids evaluation on the poles.
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.
Thanks; for some reason, I thought the grid function metadata had the appropriate spin-weight baked in.
Thanks for the time-refinement info also; I'll bear it in mind, but we have some (possibly misguided) reason for what I listed.
Beany