Generic

[znasipak]

This branch includes an implementation of Teukolsky source integration for generic geodesics in Kerr spacetime. However, before this branch is merged, we need to modify the KerrGeodesics package so that it outputs all relevant geodesic functions that are required to perform the source integration. Also, should we get rid of all of the different sources for different special cases of different orbits (e.g., spherical, eccentric vs generic)? Instead we could just have generic sources and have the ConvolveSource function handle the different limiting cases as it performs the source integration.

[barrywardell]

I've just pushed a set of updates that tidy this up and make various speed and reliability improvements. There are a couple of things I noticed/remaining:

1. Are we sure about the integration limits for the eccentric and generic cases? I thought they should be $[0, 2\pi]$ rather than $[0,\pi]$. In some cases the integrand is symmetric about $\pi$ so we can get away with it, but that's not true of all cases.
2. There is a loss of precision in the calculation of the sources near turning points. This is because we take a square root to get the four-velocity and this ends up producing a result which is half as precise as it should be.
3. The code in TeukolskySource.m could generally be tidied up now that we have the four velocity in KerrGeodesics.

The unit tests currently fail on GitHub but that is because we require a new, unreleased version of KerrGeodesics. All tests pass on my laptop.

[barrywardell]

It would be helpful to have some independent test cases to compare against. @znasipak @vskoupy @Kevin-Cunningham @MvdMeent: do you have data to compare against? Specifically, if we could check the values in the "Amplitudes" section of Tests/TeukolskyRadial.wlt that would be very helpful.

[znasipak]

@barrywardell Since switching to NIntegrate, the generic source no longer converges for any of the orbits I tried. I'm not sure if you are getting the same issue on your machine. Also, I will take a look at all of the scalar amplitudes in the tests. As for your comments above

1. As you mentioned, most terms in the source integral are symmetric under $q_r \rightarrow 2\pi - qr$ and $q\theta \rightarrow 2\pi - q_\theta$. The only terms that are not symmetric are the terms that depend on the four-velocity. So we handle those cases explicitly by summing over all orientations of $u^r$ and $u^\theta$ in the integrand itself. This way we only need to integrate from $[0, \pi]$, which reduces the number of times we have to evaluate the spheroidal harmonics and radial Teukolsky functions.
2. One way to address this issue is to add /.x_/;x==0 :> 0 for ur and u\[theta]. I'm not sure if adding this pattern matching would significantly slow down the code, though.
3. Or another way of addressing the issue may be to make use of your suggestion and incorporate the four-velocity functions introduced by KerrGeodesics

[barrywardell]

Can you give a couple of specific examples that are failing? All of the cases I tried (i.e. the ones in the unit test file) worked fine for me and were faster/more accurate/more robust than the custom integration code.

Thanks for clarifying point 1. I wondered about that after sending my message so it's great to have confirmation. Is it clear to you why when summing over orientations we still get answers that differ by more than a factor of 2 when integrating over $[0,\pi]$ rather than $[0,2\pi]$?

For point 2, there is still an issue if $u^r$ or $u^\theta$ is very close to zero so I think solution 3 might be more robust in general.

[znasipak]

@barrywardell Here is one example:

In[40]:= orbit = KerrGeoOrbit[0.9`40, 10.`40, 0.2`40, 0.9`40];

In[41]:= AbsoluteTiming[
 teuk = TeukolskyPointParticleMode[-2, 12, 2, 2, 2, orbit];]

During evaluation of In[41]:= NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in Teukolsky`ConvolveSource`Private`q\[Theta] in the region {{0,3.1416},{0,3.1416}}. NIntegrate obtained -3.9306*10^16+6.4050*10^16 I and 1.5649`4.853245542879879*^17 for the integral and error estimates.

During evaluation of In[41]:= NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in Teukolsky`ConvolveSource`Private`q\[Theta] in the region {{0,3.1416},{0,3.1416}}. NIntegrate obtained -3.8021*10^14-1.7645*10^14 I and 4.1918`4.895016300138209*^14 for the integral and error estimates.

Out[41]= {6.23733, Null}

In[42]:= teuk["Amplitudes"]

Out[42]= <|"\[ScriptCapitalI]" -> 3.5359*10^-16 - 3.1450*10^-16 I, 
 "\[ScriptCapitalH]" -> 1.9705*10^-18 + 1.7561*10^-18 I|>

In[43]:= orbitHighPrecision = 
  KerrGeoOrbit[0.9`70, 10.`70, 0.2`70, 0.9`70];

In[44]:= AbsoluteTiming[
 teukHighPrecision = 
   TeukolskyPointParticleMode[-2, 12, 2, 2, 2, orbitHighPrecision];]

During evaluation of In[44]:= NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in Teukolsky`ConvolveSource`Private`q\[Theta] in the region {{0,3.1415926535897932384626433832795},{0,3.1415926535897932384626433832795}}. NIntegrate obtained 8.2261556075019958560702924725947*10^16-4.4549902323899187980233735429900*10^17 I and 5.5639900839479195957050454883909`32.18069323465016*^6 for the integral and error estimates.

During evaluation of In[44]:= NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in Teukolsky`ConvolveSource`Private`q\[Theta] in the region {{0,3.1415926535897932384626433832795},{0,3.1415926535897932384626433832795}}. NIntegrate obtained -1.28322946290833382689879879327327*10^15-2.4894630753567939630183133621069*10^15 I and 38878.851351176722845086405747768`32.19999686537796 for the integral and error estimates.

Out[44]= {16.9003, Null}

In[45]:= teukHighPrecision["Amplitudes"]

Out[45]= <|"\[ScriptCapitalI]" -> \
-1.3072756178002441877206119929262*10^-15 + 
   2.5356115479849852859745035865568*10^-15 I, 
 "\[ScriptCapitalH]" -> 
  3.2010448049046135822100755897191*10^-18 + 
   1.73435239614377256816803960490614*10^-17 I|>

It does seem to converge to the correct answer when I throw more precision at it. However it appears that NIntegrate still fails to converge to the requested precision/accuracy. My way of interpreting this is that even the high-precision results are only accurate to 10 digits, if the error estimate warnings are accurate.

[MvdMeent]

@znasipak Are you sure that after adding together the different legs, the integrand is now a smooth pi-periodic function in both directions? If not, the Trapezoidal method will fail to converge exponentially?

[barrywardell]

@znasipak Are you sure that after adding together the different legs, the integrand is now a smooth pi-periodic function in both directions? If not, the Trapezoidal method will fail to converge exponentially?

My experiments indicated that it's not $\pi$-periodic in all cases, but it is $2\pi$ periodic.

[barrywardell]

@barrywardell Here is one example:

In[40]:= orbit = KerrGeoOrbit[0.9`40, 10.`40, 0.2`40, 0.9`40];

In[41]:= AbsoluteTiming[
 teuk = TeukolskyPointParticleMode[-2, 12, 2, 2, 2, orbit];]

During evaluation of In[41]:= NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in Teukolsky`ConvolveSource`Private`q\[Theta] in the region {{0,3.1416},{0,3.1416}}. NIntegrate obtained -3.9306*10^16+6.4050*10^16 I and 1.5649`4.853245542879879*^17 for the integral and error estimates.

During evaluation of In[41]:= NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in Teukolsky`ConvolveSource`Private`q\[Theta] in the region {{0,3.1416},{0,3.1416}}. NIntegrate obtained -3.8021*10^14-1.7645*10^14 I and 4.1918`4.895016300138209*^14 for the integral and error estimates.

Out[41]= {6.23733, Null}

In[42]:= teuk["Amplitudes"]

Out[42]= <|"\[ScriptCapitalI]" -> 3.5359*10^-16 - 3.1450*10^-16 I, 
 "\[ScriptCapitalH]" -> 1.9705*10^-18 + 1.7561*10^-18 I|>

In[43]:= orbitHighPrecision = 
  KerrGeoOrbit[0.9`70, 10.`70, 0.2`70, 0.9`70];

In[44]:= AbsoluteTiming[
 teukHighPrecision = 
   TeukolskyPointParticleMode[-2, 12, 2, 2, 2, orbitHighPrecision];]

During evaluation of In[44]:= NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in Teukolsky`ConvolveSource`Private`q\[Theta] in the region {{0,3.1415926535897932384626433832795},{0,3.1415926535897932384626433832795}}. NIntegrate obtained 8.2261556075019958560702924725947*10^16-4.4549902323899187980233735429900*10^17 I and 5.5639900839479195957050454883909`32.18069323465016*^6 for the integral and error estimates.

During evaluation of In[44]:= NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in Teukolsky`ConvolveSource`Private`q\[Theta] in the region {{0,3.1415926535897932384626433832795},{0,3.1415926535897932384626433832795}}. NIntegrate obtained -1.28322946290833382689879879327327*10^15-2.4894630753567939630183133621069*10^15 I and 38878.851351176722845086405747768`32.19999686537796 for the integral and error estimates.

Out[44]= {16.9003, Null}

In[45]:= teukHighPrecision["Amplitudes"]

Out[45]= <|"\[ScriptCapitalI]" -> \
-1.3072756178002441877206119929262*10^-15 + 
   2.5356115479849852859745035865568*10^-15 I, 
 "\[ScriptCapitalH]" -> 
  3.2010448049046135822100755897191*10^-18 + 
   1.73435239614377256816803960490614*10^-17 I|>

It does seem to converge to the correct answer when I throw more precision at it. However it appears that NIntegrate still fails to converge to the requested precision/accuracy. My way of interpreting this is that even the high-precision results are only accurate to 10 digits, if the error estimate warnings are accurate.

I also find the same thing. Is this just a particularly hard case (quite high-l, low m), or is it a case the code should be handling easily? The warning message is just a consequence of NIntegrate stopping at the default limit for MaxRecursion. We could increase this if we think it's a sensible thing to do and we have a good value to increase it to. Perhaps this is also linked to @MvdMeent's comment about exponential convergence.

[znasipak]

@znasipak Are you sure that after adding together the different legs, the integrand is now a smooth pi-periodic function in both directions? If not, the Trapezoidal method will fail to converge exponentially?

My experiments indicated that it's not π-periodic in all cases, but it is 2π periodic.

It is not $\pi$-periodic, but it is symmetric over the entire domain

[MvdMeent]

In order to get spectral convergence we need the integrand to be periodic (and smooth!) .

[znasipak]

In order to get spectral convergence we need the integrand to be periodic (and smooth!) .

Naively yes, which is why I wrote my own spectral integrator. Rewriting the integrand so that it is symmetric changes the spectral integration from taking the DFT of the integrand to the DCT. So you still get exponential convergence.

[MvdMeent]

In order to get spectral convergence we need the integrand to be periodic (and smooth!) .

Naively yes, which is why I wrote my own spectral integrator. Rewriting the integrand so that it is symmetric changes the spectral integration from taking the DFT of the integrand to the DCT. So you still get exponential convergence.

That makes sense. It sounds like the easy (quick) fix is to have the NIntegrate range be the full 2\pi. As far as I know, NIntegrate doesn't have a separate Method for integrate symmetric periodic functions. (Maybe it does?)

[znasipak]

As a check, I rewrote the source integrand and integral to be over the standard $[0,2\pi]$ domain in $qr$ and $q\theta$. You get the same results and same error messages from Mathematica. It just takes 2-4 times as long to evaluate.

[znasipak]

As a check, I rewrote the source integrand and integral to be over the standard [0,2π] domain in qr and qθ. You get the same results and same error messages from Mathematica. It just takes 2-4 times as long to evaluate.

I spoke too soon. Changing to the standard $[0, 2\pi]$ range leads to my $l=12$ case failing to converge to the right answer even with 70 digits of precision. This may be tied to the precision loss at the turning points that @barrywardell observed.

[barrywardell]

In order to get spectral convergence we need the integrand to be periodic (and smooth!) .

Naively yes, which is why I wrote my own spectral integrator. Rewriting the integrand so that it is symmetric changes the spectral integration from taking the DFT of the integrand to the DCT. So you still get exponential convergence.

That makes sense. It sounds like the easy (quick) fix is to have the NIntegrate range be the full 2\pi. As far as I know, NIntegrate doesn't have a separate Method for integrate symmetric periodic functions. (Maybe it does?)

The current version with NIntegrate uses the trapezoid rule. My understanding is that in practice the spectral integration reduces to this, is that correct? Does this hold irrespective of whether we are considering the integral as a DFT or DCT?

To check, I've just run both NIntegrate and the custom codes for an eccentric orbit with a=0.1, p=10, e=0.1, x=1. I tweaked the precision goals a little to ensure that both methods ended up evaluating the integrand the same number of times (in this case 17 evaluations). The results for the two methods agree to ~13 digits, which is consistent within the target PrecisionGoal of 12. Looking at the integrand evaluations, they are done at the same points in the two methods. See the figures below showing NIntegrate in orange and custom integration in blue. NIntegrate took ~1.1s in this case and the custom code took ~1.2s.

Note that I get a slightly different answer (differing in the second or third digit) if I integrate up $2\pi$ and divide by 2. The integrand is $2\pi$ periodic and nearly, but not quite symmetric about $\pi$. This is shown in the below plot of the integrand as a function of $q_r$ (blue) and $2\pi - q_r$ (orange).

[znasipak]

In order to get spectral convergence we need the integrand to be periodic (and smooth!) .

Naively yes, which is why I wrote my own spectral integrator. Rewriting the integrand so that it is symmetric changes the spectral integration from taking the DFT of the integrand to the DCT. So you still get exponential convergence.

That makes sense. It sounds like the easy (quick) fix is to have the NIntegrate range be the full 2\pi. As far as I know, NIntegrate doesn't have a separate Method for integrate symmetric periodic functions. (Maybe it does?)

The current version with NIntegrate uses the trapezoid rule. My understanding is that in practice the spectral integration reduces to this, is that correct? Does this hold irrespective of whether we are considering the integral as a DFT or DCT?

To check, I've just run both NIntegrate and the custom codes for an eccentric orbit with a=0.1, p=10, e=0.1, x=1. I tweaked the precision goals a little to ensure that both methods ended up evaluating the integrand the same number of times (in this case 17 evaluations). The results for the two methods agree to ~13 digits, which is consistent within the target PrecisionGoal of 12. Looking at the integrand evaluations, they are done at the same points in the two methods. See the figures below showing NIntegrate in orange and custom integration in blue. NIntegrate took ~1.1s in this case and the custom code took ~1.2s.

Note that I get a slightly different answer (differing in the second or third digit) if I integrate up 2π and divide by 2. The integrand is 2π periodic and nearly, but not quite symmetric about π. This is shown in the below plot of the integrand as a function of qr (blue) and 2π−qr (orange).

In the current code, the integrand is hard-coded to assume $q_r$ runs from $0$ to $\pi$. If you want to generalize it, change the lines

{Cnnp1p1,Cnmbarp1p1,Cmbarmbarp1p1} = TS["Cab"][r0, \[Theta]0, 1, 1];
{Cnnm1p1,Cnmbarm1p1,Cmbarmbarm1p1} = TS["Cab"][r0, \[Theta]0, -1, 1];

to

{Cnnp1p1,Cnmbarp1p1,Cmbarmbarp1p1} = TS["Cab"][r0, \[Theta]0, Sign[Sin[qr]], 0];
{Cnnm1p1,Cnmbarm1p1,Cmbarmbarm1p1} = TS["Cab"][r0, \[Theta]0, -Sign[Sin[qr]], 0];

[barrywardell]

Ah yes, that does restore symmetry, and produces the same answer as the more efficient half-interval version.

[barrywardell]

Getting back to @znasipak's test case, there seem to be two issues:

1. NIntegrate is complaining about not converging to the specified precision.
2. It's a difficult case susceptible to roundoff cancellations.

I haven't been able to get either integration (NIntegrate or custom) to produce the right answer when working at 40-digit precision; this is probably due to point 2. Switching to 70 digits, both methods produce consistent answers. I find that NIntegrate is ~5 times faster in this case, but that is mostly due to it only doing ~1/4 of the number of integrand evaluations before it decides it is not converging further.

[MvdMeent]

It would be helpful to have some independent test cases to compare against. @znasipak @vskoupy @Kevin-Cunningham @MvdMeent: do you have data to compare against? Specifically, if we could check the values in the "Amplitudes" section of Tests/TeukolskyRadial.wlt that would be very helpful.

I've checked the the s=-2 amplitudes against my code and they agree.

[barrywardell]

The loss of precision at the turning points is now fixed for the s=-2 generic case. The same thing still needs to be done for the other cases.

[znasipak]

It would be helpful to have some independent test cases to compare against. @znasipak @vskoupy @Kevin-Cunningham @MvdMeent: do you have data to compare against? Specifically, if we could check the values in the "Amplitudes" section of Tests/TeukolskyRadial.wlt that would be very helpful.

The scalar amplitudes agree (to within 12+ digits for my C++ code results). The only issue is that for eccentric orbits the In and Up values are switched. I'll create an issue for this