astro-group-bristol / Gradus.jl

Extensible spacetime agnostic general relativistic ray-tracing (GRRT).
https://astro-group-bristol.github.io/Gradus.jl/dev/
GNU General Public License v3.0
18 stars 2 forks source link

Lineprofiles: better handling of initial guess in root finder #72

Closed phajy closed 1 year ago

phajy commented 1 year ago

I was going to post this as a follow-up in Precision tolerances but have decided to make it a "new" issue instead.

I've encountered a (possibly?) precision-related error: Roots.ConvergenceFailed("Algorithm failed to converge"). Here's an example to investigate (based on this example).

a = 0.0
d = GeometricThinDisc(0.0, 400.0, π / 2)
u = @SVector [0.0, 1000.0, deg2rad(60), 0.0]
m = KerrMetric(1.0, a)

# maximal integration radius
maxrₑ = 200.0

# emissivity function
ε(r) = r^(-3)

# g grid to do flux integration over
gs = range(0.0, 1.2, 500)
_, flux = lineprofile(gs, ε, m, u, d, maxrₑ = maxrₑ)

I've briefly experimented with some of the hidden parameters, e.g., increasing the number of g* points or reducing the tolerance but with no success. This differs from the example by having a much larger outer radius where $g\text{max} - g\text{min}$ will be very small and close to zero.

fjebaker commented 1 year ago

I am able to reproduce. These sorts of precision errors in the root finder normally mean the extrema of the bracketing interval which sets up the precision tracer missed the disc. This often happens when roughly $r\text{e} \geq \frac{1}{2} r\text{d}$.

A simple fix is to set, e.g.

d = GeometricThinDisc(0.0, 800.0, π / 2)

But I'll leave this issue open with the addition that

For the latter, the offending line is

https://github.com/astro-group-bristol/Gradus.jl/blob/ec031bb31a9458a97816859b5b52257cb5102e40/src/transfer-functions/cunningham-transfer-functions.jl#L175

which could probably be safely set to d.outer_radius + Δ, where Δ is just a safety offset. In the past this had to scale to help the root finder converge, but we're using a better default algorithm now, so we might be able to get away with just using the same maximum offset for each emission radius.

fjebaker commented 1 year ago

Fixed in #107