[Users] using WeylScal4 and Multipole for WF output -- scheduling issues

[Ian Hinder]

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.

-- 
Ian Hinder
ian.hinder at aei.mpg.de