astro-group-bristol / gradus.jl Goto Github PK

Use a quad-tree to make resolutions of any 2d plot better near distorted edges better, e.g. for exclusions zones in naked singularities, or disc edges or close to the event horizon in rendering.

Would work similarly to adaptive refinement, and could maybe even vendor a library for it?

Tracking allocations and analyzing with Coverage.jl shows that we allocate enormously during problem setup:

 CoverageTools.MallocInfo(0, "dev/Gradus/src/metrics/boyer-lindquist-fo.jl.101246.mem", 335)
 CoverageTools.MallocInfo(736, "dev/Gradus/src/Rendering/render.jl.101246.mem", 80)
 CoverageTools.MallocInfo(736, "dev/Gradus/src/Rendering/render.jl.101246.mem", 81)
 CoverageTools.MallocInfo(6528, "dev/Gradus/src/Rendering/render.jl.101246.mem", 77)
 CoverageTools.MallocInfo(7360, "dev/Gradus/src/GeodesicTracer/constraints.jl.101246.mem", 29)

The matrix dispatch already has inplace modifications, but for e.g. vectors of static vectors:

There is no need to use map here, we should be able to re-use the same container.

The top row of an image calculated with rendergeodesics is actually the bottom row.

Keyword arguments are passed to the tracer in trace_equitorial_circular_orbit, since upper_rate and lower_rate are not seperately captured / kwargs passed to both the trace function and solve function.

This is probably related to reusing the integrator state in each thread, but e.g. Feagin10 or Vern6 fail dramatically when used in the context of e.g. lineprofiles, and I think this is the reason.

The disc intersection algorithm is poor, and may return very different results for slightly different initial parameters (hinted to by the tearing effects at the edge of the disc for large step sizes).

This really needs to be fixed, as it is impacting tolerances for all subsequent calculations.

Using this issue to track progress on clamping down why the transfer function integration and the regular binning method for calculating the line profiles disagree.

Crux of the problem:

Code to reproduce

using Gradus
using Plots

function calculate_line_profile_binned(m, x, d, bins; maxrₑ = 500.0, minrₑ = Gradus.isco(m) + 1e-2)
    plane = PolarPlane(GeometricGrid(); Nr = 1000, Nθ = 1000, r_max = 2 * maxrₑ)
    _, f = lineprofile(
        m, 
        x, 
        d, 
        algorithm = BinnedLineProfile(), 
        verbose = true,
        bins = bins,
        maxrₑ = maxrₑ,
        minrₑ = minrₑ,
        plane = plane
    )
    f
end

function calculate_line_profile(m, x, d, bins; maxrₑ = 500.0, minrₑ = Gradus.isco(m) + 1e-2) 
    _, f = lineprofile(
        m, 
        x, 
        d, 
        algorithm = CunninghamLineProfile(), 
        verbose = true,
        bins = bins,
        maxrₑ = maxrₑ,
        minrₑ = minrₑ,
        numrₑ = 200,
        β₀ = 2,
        Nr = 3000,
    )
    f
end

bins = collect(range(0.1, 1.4, 200))
m = KerrMetric(1.0, 0.998)
x = SVector(0.0, 1000.0, deg2rad(45), 0.0)

ssd = ShakuraSunyaev(m, eddington_ratio = 0.3)

tf30 = @time calculate_line_profile(m, x, ssd, bins)
bn30 = @time calculate_line_profile_binned(m, x, ssd, bins)

plot(bins, tf30, label = "transfer functions")
plot!(bins, bn30, label = "binned")

There already is a function in GradusBase defined for mapping vectors to local sky angles

but this is only implmented for the Kerr case at the moment. We either would need a good way of calculating this in other metrics, else find a reliable way to do parallel transport.

Create a vector wrapper type that knows which coordinates it is defined in; essentially an SVector with some additional parameter that encodes if it's a spherical or cartesian vector.

This would allow much more expressive coordinate transformation functions, and fewer mistakes made by users.

Discussed in #73

Continuing discussion in this issue.

Looking at the profiler, we might be able to squeeze performance by evaluating the metric_components function along with the Jacobian to save a little bit of time. I think there's a way to do this with ForwardDiff.jl, but opening an issue so I remember to look at this later:

A test runner for PRs should be opened that ensures the examples still pass, as they occasionally break without breaking smoke tests (which also implies the smoke test coverage is poor), e.g. #76

Makes sense doesn't it. They are the same thing, and in cylindrical coordinates, the (projected) radius is denoted $\rho$.

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.

Modify cross_section and related functions to receive only $\rho$ instead of the full four vector. There could be additional structures, like AsymmetricThickDisc that permit more complex structures, such as thick discs that change in the $\phi$ coordinate.

With regard to GC, I noticed that lagtransfer can spend a lot of time no garbage collection. 3533.139993 seconds (292.44 M allocations: 87.741 GiB, 83.50% gc time, 0.02% compilation time).

-- @phajy in #30

This is an insane amount of time lost to GC and needs to urgently be addressed.

This has a number of things to it

• Update doc strings for all functions and structures
• Write introduction and about pages
• Write examples pages
• Write method detail pages
• Compile full API docs

This issue is carried over from when we used SageMath to generate the symbolic geodesic equations:

• fjebaker/ComputedGeodesicEquations.jl#18

JP metric has the sign of the spin reversed, e.g. a=1 corresponds to a=-1 for Kerr.

The issue was not the sign convention as mentioned in the issue above, but is related to using the wrong sign in the constraint equation.

I am yet to check if the AD methods also suffer from this, but leaving the issue open until I've had a chance to confirm.

If the outer radius on the image plane is not a few factors of maxrₑ, which may occur if a user is specifying the plane, then there is no warning that the sampled region may not actually be the desired region.

The image_plane function should be reworked, so that the planes themselves store parameters of the plane, but the size of it is scaled internally, so that this mistmatch cannot occur.

This should also still be useable in the rendergeodesic functions however, so some care in the design needs to be taken.

As part of our tests, we should also keep a benchmark suite to track the performance of our code, and make sure new changes don't accidentally cause our integration to slow down.

This should have --math-mode=ieee and a variety integration tolerances and image sizes to ensure our methods scale well, and don't git the GC too hard (as they have in the past).

The ThinDisc has the OpticallyThin trait.

I doubt that the relativistic beaming is being properly accounted for, since we do not make a Lorentz transform before the tetrad transformation.

It would be better to get the transfer functions to sample evenly in $g^\ast$ instead of in $\theta$.

The proposed algorithm change is

• estimate gmin gmax coarsely.
• for $g^\ast$ in set of target $g^\ast$, root find for radius, then root find over $\theta$ for the corresponding $g^\ast$.

An issue that might crop up here are those transfer functions where the inferred extrema are worse than the area extrema. In these cases we normally see an error of 1e-4 or so, but that could still be quite damaging to the calculation. Perhaps, however, we can spend more time in the extrema optimizer at the trade off of needing fewer $g^\ast$ than theta to sample the transfer functions fully.

• position to be x, velocity to be v, disc velocities to be u.
• standard convention should be to have observer at 10'000 $r_\text{g}$, not 1000.

At the moment the user has to do some pretty non-trivial arithmetic to work out what the corresponding impact parameter axes are when using the *rendergeodesics methods. I propose instead that a new API works like this:

a, b, img = rendergeodesics(...)

where a and b are the $\alpha$ and $\beta$ impact parameters respectively.

Could maybe make it easier to switch between pinhole and photographic plane observers, since at the moment we're hardcoding pinhole observer.

To facilitate a nicer analysis API, I suggest we move to wrap the SciMLBase types in our own struct, so that we can dispatch methods more targetedly. This would allow us to add plotting recipes unique to our package, allow overloads with the statistics community's fit functions, and generally let us be more expressive.

We could then also track which coordinates a given solution is in (Cartesian vs spherical), and provide mappings between them.

Takes about two minutes to pre-compile. This should really be reduced, there is no reason for this to be so long.

See comments in:

• #28

MacOS seems to produce wildly different values in the tests in the CI. Need to find out why, and then add other OSs.

There is a relatively small chance that this happens:

    nested task error: Something went wrong. Integrator stepped past tstops but the algorithm was dtchangeable. Please report this error.
    Stacktrace:
     [1] error(s::String)
       @ Base ./error.jl:35
     [2] handle_tstop! ...

Given that it is non-deterministic, it is almost certainly related to some threading issue and a race condition being hit.

Turns out this is all related to weird overflows with the progress meter. Commenting out any and all progress meter calls fixes everything, so we need some sort of thread safety here / it is likely that copying the meter state per threads is the culprit.

I'm trying to calculate the ISCO for a black hole spinning at (almost) the limit allowed by the spacetime. Perhaps there is some very small tolerance issue with the integration? It is going to be numerically difficult to figure out both 1) what is the maximum allowed value of $a$, and 2) what is the corresponding ISCO.

Here's an example

m = JohannsenPsaltisMetric(M = 1.0, a = 0.8168816881688169, ϵ3 = 0.8)
is_naked_singularity(m)
Gradus.isco(m)

says that the singularity is not naked but when calculating the ISCO gives an error

ERROR: ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.

I've labelled this "bug" but perhaps it isn't really. If I make $a$ just ever so slightly smaller (e.g., by $\sim 10^{-10}$ [!!]) it works. So maybe my $a$ is not really valid and the singularity is naked for that value (might be hard to determine because at the limit the singularity is "only just" naked).

Not urgent!

Self descriptive

For deviations from 90 degress in theta, the beta impact parameter gets stretched, distorting the aspect ratio of the disc. Most noticeable for <60 degrees.

apply is returning a flat array instead of matching the behaviour of rendergeodesics with a point function.

There's several arbitrarily set precisions in the tests, and more generally some unusually high discrepencies in the values calculated between different methods (see brief mention in this blog post).

We should have a suite of tests that

• ensures the precisions of all methods are known
• attempts to minimize any extraneous errors
• keeps the different integration methods consistent with eachother

I think as an aim, we want our methods to agree within 1e-8 or 1e-9 for relative tolerances. In the literature this isn't really discussed other than occasionally ensuring

$$ g_{\sigma\nu} \dot{u}^\sigma \dot{u}^\nu + \mu^2 < \varepsilon \approx 10^{-8} $$

and even that has a big "-ish" stamped onto the end, as some authors don't even really seem to check that, and others just check whether it agrees with some other simulation.

I think we can do better and quantify the errors more.

Update the getting started guide and docstrings to reflect the definition of impact parameters

$$ \alpha := \frac{ p_{(\phi)} }{ p_{(t)} }, \quad \beta := \frac{ p_{(\theta)} }{ p_{(t)} } $$

and how they are mapped from an (stationary) LNRF.

The default value of extrapolate in DataInterpolations has been changed to false. See SciML/DataInterpolations.jl#194.

This can, in some cases, break Gradus.

Adding ;extrapolate=true) to all of the DataInterpolations.LinearInterpolation calls in the following lines seems to fix the problem.

I would guess that extrapolations near the edge of the domain are expected. Not sure if this needs to be looked at a bit more?

The project is a little over a year old now and the code bloat is apparent.

Need to go through each function with a fine toothed comb and see what can be made leaner, or even just removed entirely. Also there are a number of bad habits that should be corrected, such as

• remove always using where {T} unless actually needed
• anonymous lambda functions are difficult to debug, give them names
• function names and extraction API, e.g. getgeodesicpoint vs just process, as that is what it's really doing.
• #16

Also worth looking into

• #13
• #8 : not worth it as it currently stands

before getting too carried away with more features.

Currently, there are some implemenations that are

metric_components(m::AbstractMetric, x::SVector{4})

and others that are

metric_components(m::AbstractMetric, xt::SVector{2}) # only r and theta components
metric_components(m::AbstractMetric, xt::NTuple{2}) # same as above

The return type is also all over the place.

This is needlessly confusing. Suggest we differentiate metric_components and metric_elements, the prior is used as the implementation, and metric_elements is used as an accessor method?

This also seems complex though...

Instead of pre-generating the whole u and v array, there should be a dispatch of tracegeodesic where both position and velocity are functions, such that these can be generated on the fly to avoid un-needed allocations.

As described in the comments of

• #71

when the observer angle is close to 85, the root finder misses the disc underneath the equitorial plane.

Code to reproduce (v0.2.3):

using Gradus, StaticArrays, Plots
m = KerrMetric(M=1.0, a=0.998)
d = GeometricThinDisc(0.0, 3000.0, π/2)

angle = 85

u = @SVector [0.0, 1000.0, deg2rad(angle), 0.0]

# errors here
ctf = cunningham_transfer_function(m, u, d, 4.0)

mask = @. (ctf.g✶ > 0.001) & (ctf.g✶ < 0.999) 
plot(ctf.g✶[mask], ctf.f[mask], label="θᵢ=$angle", xlabel="g✶", ylabel="f") 

Setting e.g. offset_max = 15.0 seems to fix it, but produces visible errors in the resulting transfer function:

This was previously discussed in

• fjebaker/GeodesicRendering.jl#9

but updating the state of the problem.

The render caches are quite enormous, much larger than they really should be at 280MB for a 1250x1750 render. I mentioned writing a custom serialiser / de-serialiser to this end for JLD2, which is still an option.

Perhaps we also look into StructArrays.jl for this? May help with storing the geodesic point types, and performance, since we commonly access similar fields (e.g., position).

Opening an issue so I remember to address this later: when using the (forthcoming) LNRF calculator, the tests for the Dilaton-Axion metric fail, since

4×4 Matrix{Float64}:
 1.0          0.0  0.0   0.0
 0.0          1.0  0.0   0.0
 0.0          0.0  1.0   0.0
 2.22045e-16  0.0  0.0  -1.0

is not Minkowski. We have the correct signature, but it seems to be in the wrong order?

It's been a few weeks since the last release, and much has changed in the main branch. This issue is to track what still needs to be done before 0.4.17 0.5.0 release:

• Documentation examples need to be updated and recomputed (#168 )
• #154 needs to be documented
• #150 needs to be documented
• #148 needs to be documented (debug toolkit)
• #157 needs to be documented
• #141 needs to be documented
• #158 needs to be documented (breaking change)
• #109
• #156
• #161 needs to be documented
• #160
• #110
• #166

From #143

The seam in this figure corresponds to the bin edge $g^\ast \in [0, h]$ and $g^\ast \in [1-h, 1]$. I tried a few different schemes to resolve this, but I think it's a combination of getting the timing wrong, and also the flux in the bin edges being too low with that Dauser approximation?

The gtol parameter is a mess. Trying to do something "sophisticated" like in #172 has led to terrible impulse response functions:

Using the old linear gtol * r scaling is already much much better:

Looking at the 2d integrated heatmap:

"Tuning" (read: randomizing) the magic constants in the _gtol_error function gives this

This is most likely due to Jacobians going all over the place at the boundaries of _gtol_error, which means we either need any old $C^\infty$ function, or just go back to the linear scaling.

I'm really not happy with how the transfer functions are currently defined and used. Since the integration domain is over $g$ when calculating the line profiles anyway, the conversion to $g^\ast$ is really inconvenient, as it introduces a number of numerical issues when calculating the transfer functions themselves. Consider the definition

$$ f := \frac{g}{\pi r} \sqrt{g^\ast (1 - g^\ast)} \left\lvert \frac{\partial(\alpha \beta}{\partial(r, g^\ast)} \right\rvert, $$

which is then used to find the flux through

$$ \text{d}F = g^2 I \frac{\pi r}{ \sqrt{g^\ast (1 - g^\ast)} } f \text{d}r \text{d} g^\ast. $$

This is itself used in the integral

$$ F = \int_{r=r_0}^{r=r_\text{max}}\int_{g^\ast=0}^{g^\ast=1} \text{d}F. $$

It's pretty obvious that most of the terms cancel. Care when integrating comes from the regions where $g^\ast \rightarrow { 0,1 }$, but that is because by the transfer functions go to zero at those points, but $\text{d}F \rightarrow \infty$ numerically, even though there is an analytic value. The prescription is then to have a buffer $h &gt; 0$, where special edge integrations are performed instead.

The transfer functions could far more conveniently be defined

$$ f := g\left\lvert \frac{\partial(\alpha \beta}{\partial(r, g)} \right\rvert, $$

such that

$$ \text{d}F = g^2 I f \text{d}r \text{d} g, $$

and the integration is then performed directly over $g$ instead of $g^\ast$, which also then elides the need for a $g \rightarrow g^\ast$ Jacobian.

Possible gripes can easily be fixed in this formalism too:

• With $g^\ast$ we can integrate over the full energy. So too with $g$, since we need to rebin anyway when extracting the line profiles, so we just integrate over the limits of the bins.
• How would one know when to use which transfer function? Easy, simply set $f=0$ everywhere except $g_\text{min} \leq g \leq g_\text{max}$. Infact, one can even simplify the total integral this way too, by adding the transfer functions all together for different radi, and then integrating over $g$.

It seems the transfer functions could be so much easier to use when defined this way.

The equitorial disc velocities used to calculate the redshift are hard coded for the Kerr metric. This should be seperated out, and then allowed to be made generic.

Additionally, should try and build in some generic analytic solutions for e.g. all static, axis-symmetric metrics.

Applies to the redshift calculations in the plunging region. See discussion in:

• astro-group-bristol/AccretionFormulae.jl#22

parity of the radial component of velocity seems to be opposite under time reversal

We need a justification for this.

This has been a point of contention for a while. The current state of DiffEqGPU.jl does not support SecondOrderODEProblems out of the box, and a patch I wrote for DiffEqGPU's problematic parse reveals a more fundemental problem:

• The full integration batch terminates when a single geodesic terminates. If a single geodesic falls into the black hole or intersects with the disc, the whole simultaneous solve is terminated.

We also need to evaluate to what degree we want to provide GPU support for the integration problems, as it may reduce us to Float32 operations, and be an awful lot of work for marginal performance gain.

A single trace is a relatively cheap problem, in the micro-second solve time, so it should be entirely feasible and logical to chuck this on the GPU -- the question is just how.

As a note, the above issue with SecondOrderODEProblems can also be solved by converting our problem to a first order one and manually updating the position vector.