On 22 Oct 2010, at 15:32, Bernard Kelly wrote:
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.
Scheduling in Cactus, especially with mesh refinement, is very complicated. The schedule for WeylScal4 has been carefully chosen to work (by me, Roland and Erik) - if you change it, you're on your own :)
Why did that change?
The difference is related to prolongation of mesh refinement boundaries. Psi4 is computed using spatial derivatives, and hence needs the boundaries of the grid to be filled in via prolongation from the next coarser grid. Functions scheduled in ANALYSIS cannot have time interpolation, due to the time at which ANALYSIS is scheduled relative to the Berger-Oliger time stepping. In order that all points in the grid have sensible values, analysis quantities which take spatial derivatives should be scheduled in EVOL so that time interpolation onto their refinement boundaries is possible. It is also necessary to recompute these quantities after regridding and in several other places. To make this convenient, a new schedule group, MoL_PseudoEvolution, was introduced in EVOL in which such quantities should be computed. Variables set in this group should have three timelevels, just like evolved variables, since time interpolation requires this.
Erik: is this description accurate?
Are variables in MoL_PseudoEvolution available for ANALYSIS automatically (i.e. will they retain their stored values)?
This depends on how you schedule storage for the variables. If you schedule storage for the variables unconditionally, then yes. In fact, since output of the variables will happen after ANALYSIS, the storage had better be permanent otherwise you won't be able to output them.
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.
We used to not have mesh refinement :)
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 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.
WeylScal4_Calculate is no longer scheduled "as calc_np" as it used to be, so there is no dependency information between it and Multipole. In the unmodified version, the dependency was automatically enforced because WeylScal4 was scheduled in MoL_PseudoEvolution and Multipole was scheduled in ANALYSIS.
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?
When Multipole calls the interpolator to find the values of Psi4 on the extraction spheres, Carpet uses the finest grid that the point exists on. This could mean the last points in a buffer zone, and if you haven't filled these correctly with prolongation, these will be poisoned (or undefined, if you're not using poisoning). So if your extraction spheres intersect your refinement boundaries (very difficult to avoid in 3D), you might pick up poison from them.
In my opinion, the interpolator should not give points from buffer zones, and instead use the coarser grid. I think this might be the case in the new (Mercurial) version of Carpet. Erik: is that the case?
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.
Right, or compute sin(th) Psi4, taking the appropriate limit in the neighbourhood of the axis. But I was more interested in computing the gravitational waves themselves, not the mode decomposition. At Scri+, Psi4 is supposed to represent the waves, a physical observable, right?
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.
No - the spin weight is not stored in the grid function. But that's a good idea!