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

[Ian Hinder]

On 22 Oct 2010, at 17:09, Bernard Kelly wrote:

> 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).

Sounds like the symmetry boundaries are not being applied.  Aha!  You  
need to set the parameters

	WeylScal4::Psi4r_group_bound = "flat"
	WeylScal4::Psi4i_group_bound = "flat"

Have a look at the test suite parameter files for WeylScal4.  This is  
not necessary in the current development branch, as this is done  
automatically by Kranc now.

It occurs to me that we could do with an example for the use of  
WeylScal4 and Multipole, e.g. by using WeylScal4 in one of the  
McLachlan QC0 example parameter files.



>
> 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
> -----------------------------------------------
>

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