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

[Bernard Kelly]

Replying to myself on this ...

I reverted the WeylScal4 schedule.ccl to its repository version didn't 
fix things -- still NaNs almost all the time.

I realised I wasn't accurate when I said I saw NaNs. What I actually see 
is that at t=0, components extracted
at r = 45, 50, 55, and SOMETIMES 60 have finite values. At all later 
times, all modes are NaN. All r=65, 70 modes
are always NaN, even at t=0.

I've tried your "offset" parameter, but that doesn't seem to help, I'm 
afraid (actually, it's on by default anyway). I've
just enabled Multipole's "out_1d_every" option, and in fact the data has 
nan at the -equator- (that is, at points
theta = 1.530519 & 1.611073 -- those closest to the x-y plane).

Bernard
> 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:
>>
>> 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?
> 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
>


-- 
-----------------------------------------------
Bernard J. Kelly

NASA Goddard Space Flight Center, Code 663
8800 Greenbelt Road
Greenbelt, MD 20771, U.S.A.

phone: +1 (301) 286-7243
-----------------------------------------------