From 66e3b90b331dcb25dfe2db6f7e3c3dfe2ce74558 Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Fri, 22 Apr 2022 19:10:53 +0100 Subject: [PATCH 1/8] added Optim --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ffb140ea..655b2b70 100644 --- a/Project.toml +++ b/Project.toml @@ -11,10 +11,11 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GeometricalPredicates = "fd0ad045-b25c-564e-8f9c-8ef5c5f21267" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d" +Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" From 9160b86859a403d6317e107003ec6f38795be40a Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Fri, 22 Apr 2022 19:11:09 +0100 Subject: [PATCH 2/8] automated circular orbits --- src/AccretionFormulae/AccretionFormulae.jl | 7 +- src/AccretionFormulae/orbit-discovery.jl | 108 +++++++++++++++++++++ 2 files changed, 114 insertions(+), 1 deletion(-) create mode 100644 src/AccretionFormulae/orbit-discovery.jl diff --git a/src/AccretionFormulae/AccretionFormulae.jl b/src/AccretionFormulae/AccretionFormulae.jl index 1a3cf128..576252d9 100644 --- a/src/AccretionFormulae/AccretionFormulae.jl +++ b/src/AccretionFormulae/AccretionFormulae.jl @@ -6,10 +6,15 @@ import ..GradusBase: AbstractMetricParams, metric using ..FirstOrderMethods: FirstOrderGeodesicPoint using ..Rendering: PointFunction +using Optim: optimize, minimizer, GoldenSection, Brent +using DocStringExtensions +using StaticArrays + include("redshift.jl") +include("orbit-discovery.jl") const redshift = PointFunction(_redshift_guard) -export redshift +export solve_equitorial_circular_orbit, trace_equitorial_circular_orbit, redshift end # module diff --git a/src/AccretionFormulae/orbit-discovery.jl b/src/AccretionFormulae/orbit-discovery.jl new file mode 100644 index 00000000..dd874d65 --- /dev/null +++ b/src/AccretionFormulae/orbit-discovery.jl @@ -0,0 +1,108 @@ + +""" + $(TYPEDSIGNATURES) + +Quality of stability measure, which has a minima for circular orbits. Effectively +a sum of the normalised residuals. +""" +Qs(rs) = sqrt(sum((rs ./ rs[1] .- 1.0) .^ 2) / length(rs)) + +function trace_single_orbit(m, r, vϕ; max_time=300.0, μ=1.0, tracer_args...) + # fixed in equitorial plane + u = @SVector[0.0, r, deg2rad(90.0), 0.0] + v = @SVector[0.0, 0.0, 0.0, vϕ] + Gradus.tracegeodesics( + m, + u, + v, + (0.0, max_time); + μ = μ, + tracer_args... + ) +end + +function measure_stability(m::AbstractMetricParams{T}, r, vϕ; tracer_args...) where {T} + sol = trace_single_orbit(m, r, vϕ; tracer_args...) + rs = selectdim(sol, 1, 6) + Qs(rs) +end + +function __solve_equitorial_circular_orbit( + m::AbstractMetricParams{T}, + r, + optimizer, + lower_bound, + upper_bound; + tracer_args... +) where {T} + res = optimize( + vϕ -> measure_stability( + m, + r, + vϕ + ; + tracer_args... + ), + lower_bound, + upper_bound, + optimizer + ) + minimizer(res) +end + +function solve_equitorial_circular_orbit( + m::AbstractMetricParams{T}, + r::Number; + lower_bound = 0.0, + upper_bound = 1.0, + optimizer = GoldenSection(), + tracer_args... +) where {T} + __solve_equitorial_circular_orbit( + m, r, optimizer, lower_bound, upper_bound + ) +end + +function sliding_window( + func, N, lower_bound, upper_bound, + lower_rate, upper_rate +) + map(1:N) do i + res = func((i, lower_bound, upper_bound)) + lower_bound = res * lower_rate + upper_bound = res * upper_rate + res + end +end + + +function solve_equitorial_circular_orbit( + m::AbstractMetricParams{T}, + r_range::Union{AbstractRange,AbstractArray}; + lower_bound = 0.0, upper_bound = 1.0, + lower_rate = 0.98, upper_rate = 1.5, + kwargs... +) where {T} + r_range_reverse = sort(r_range) |> reverse + candidate_vϕ = sliding_window( + length(r_range_reverse), + lower_bound, upper_bound, + lower_rate, upper_rate + ) do (i, lower_bound, upper_bound) + r = r_range_reverse[i] + solve_equitorial_circular_orbit(m, r, lower_bound = lower_bound, upper_bound = upper_bound) + end + reverse!(candidate_vϕ) +end + +function trace_equitorial_circular_orbit( + m::AbstractMetricParams{T}, rs; kwargs... +) where {T} + map(zip( + rs, + solve_equitorial_circular_orbit(m, rs; kwargs...) + ) + ) do (r,vϕ) + trace_single_orbit(m, r, vϕ; kwargs...) + end +end \ No newline at end of file From 8f92261940170b0b0b93500a0d862137ad535b34 Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Fri, 22 Apr 2022 19:11:25 +0100 Subject: [PATCH 3/8] function signatures for other special radii --- src/special-radii.jl | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/src/special-radii.jl b/src/special-radii.jl index c3e8eb07..08eaebfc 100644 --- a/src/special-radii.jl +++ b/src/special-radii.jl @@ -1,6 +1,27 @@ """ -$(TYPEDSIGNATURES) + $(TYPEDSIGNATURES) Innermost stable circular orbit. """ isco(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).") + +""" + $(TYPEDSIGNATURES) + +Photon orbit. +""" +r_ph(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).") + +""" + $(TYPEDSIGNATURES) + +Marginally bound orbit. +""" +r_mb(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).") + +""" + $(TYPEDSIGNATURES) + +Event horizon. +""" +r_s(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).") From 517504fa65b7fd6ac3343b91163e606dbd47ff87 Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Fri, 22 Apr 2022 19:11:43 +0100 Subject: [PATCH 4/8] update exports --- src/Gradus.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Gradus.jl b/src/Gradus.jl index 778bca6d..80b27a4d 100644 --- a/src/Gradus.jl +++ b/src/Gradus.jl @@ -73,6 +73,6 @@ include("special-radii.jl") include("AccretionFormulae/AccretionFormulae.jl") using .AccretionFormulae -export redshift +export solve_equitorial_circular_orbit, trace_equitorial_circular_orbit end # module From 712f949f1e9dd55e65aff5ae61ec2ccf7266d3ef Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Fri, 22 Apr 2022 19:11:58 +0100 Subject: [PATCH 5/8] set default tolerance to 1e-9 --- src/GeodesicTracer/GeodesicTracer.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/GeodesicTracer/GeodesicTracer.jl b/src/GeodesicTracer/GeodesicTracer.jl index 4b6bfd52..6792decc 100644 --- a/src/GeodesicTracer/GeodesicTracer.jl +++ b/src/GeodesicTracer/GeodesicTracer.jl @@ -61,8 +61,8 @@ function tracegeodesics( time_domain, solver; callback = nothing, - abstol = 1e-8, - reltol = 1e-8, + abstol = 1e-9, + reltol = 1e-9, solver_opts..., ) end From b5c112b0fcaf78d3bfd0af24f217af8f03eee0d3 Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Fri, 22 Apr 2022 19:14:49 +0100 Subject: [PATCH 6/8] julia formatting --- src/AccretionFormulae/orbit-discovery.jl | 75 ++++++++++-------------- 1 file changed, 32 insertions(+), 43 deletions(-) diff --git a/src/AccretionFormulae/orbit-discovery.jl b/src/AccretionFormulae/orbit-discovery.jl index dd874d65..a9721e6a 100644 --- a/src/AccretionFormulae/orbit-discovery.jl +++ b/src/AccretionFormulae/orbit-discovery.jl @@ -7,18 +7,11 @@ a sum of the normalised residuals. """ Qs(rs) = sqrt(sum((rs ./ rs[1] .- 1.0) .^ 2) / length(rs)) -function trace_single_orbit(m, r, vϕ; max_time=300.0, μ=1.0, tracer_args...) +function trace_single_orbit(m, r, vϕ; max_time = 300.0, μ = 1.0, tracer_args...) # fixed in equitorial plane - u = @SVector[0.0, r, deg2rad(90.0), 0.0] - v = @SVector[0.0, 0.0, 0.0, vϕ] - Gradus.tracegeodesics( - m, - u, - v, - (0.0, max_time); - μ = μ, - tracer_args... - ) + u = @SVector [0.0, r, deg2rad(90.0), 0.0] + v = @SVector [0.0, 0.0, 0.0, vϕ] + Gradus.tracegeodesics(m, u, v, (0.0, max_time); μ = μ, tracer_args...) end function measure_stability(m::AbstractMetricParams{T}, r, vϕ; tracer_args...) where {T} @@ -33,19 +26,13 @@ function __solve_equitorial_circular_orbit( optimizer, lower_bound, upper_bound; - tracer_args... + tracer_args..., ) where {T} res = optimize( - vϕ -> measure_stability( - m, - r, - vϕ - ; - tracer_args... - ), + vϕ -> measure_stability(m, r, vϕ; tracer_args...), lower_bound, upper_bound, - optimizer + optimizer, ) minimizer(res) end @@ -53,20 +40,15 @@ end function solve_equitorial_circular_orbit( m::AbstractMetricParams{T}, r::Number; - lower_bound = 0.0, + lower_bound = 0.0, upper_bound = 1.0, optimizer = GoldenSection(), - tracer_args... + tracer_args..., ) where {T} - __solve_equitorial_circular_orbit( - m, r, optimizer, lower_bound, upper_bound - ) + __solve_equitorial_circular_orbit(m, r, optimizer, lower_bound, upper_bound) end -function sliding_window( - func, N, lower_bound, upper_bound, - lower_rate, upper_rate -) +function sliding_window(func, N, lower_bound, upper_bound, lower_rate, upper_rate) map(1:N) do i res = func((i, lower_bound, upper_bound)) lower_bound = res * lower_rate @@ -79,30 +61,37 @@ end function solve_equitorial_circular_orbit( m::AbstractMetricParams{T}, r_range::Union{AbstractRange,AbstractArray}; - lower_bound = 0.0, upper_bound = 1.0, - lower_rate = 0.98, upper_rate = 1.5, - kwargs... + lower_bound = 0.0, + upper_bound = 1.0, + lower_rate = 0.98, + upper_rate = 1.5, + kwargs..., ) where {T} r_range_reverse = sort(r_range) |> reverse candidate_vϕ = sliding_window( - length(r_range_reverse), - lower_bound, upper_bound, - lower_rate, upper_rate + length(r_range_reverse), + lower_bound, + upper_bound, + lower_rate, + upper_rate, ) do (i, lower_bound, upper_bound) r = r_range_reverse[i] - solve_equitorial_circular_orbit(m, r, lower_bound = lower_bound, upper_bound = upper_bound) + solve_equitorial_circular_orbit( + m, + r, + lower_bound = lower_bound, + upper_bound = upper_bound, + ) end reverse!(candidate_vϕ) end function trace_equitorial_circular_orbit( - m::AbstractMetricParams{T}, rs; kwargs... + m::AbstractMetricParams{T}, + rs; + kwargs..., ) where {T} - map(zip( - rs, - solve_equitorial_circular_orbit(m, rs; kwargs...) - ) - ) do (r,vϕ) + map(zip(rs, solve_equitorial_circular_orbit(m, rs; kwargs...))) do (r, vϕ) trace_single_orbit(m, r, vϕ; kwargs...) end -end \ No newline at end of file +end From 03d3e87801e240a9308a43462ee91699ae84f63d Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Fri, 22 Apr 2022 19:25:38 +0100 Subject: [PATCH 7/8] circular orbit smoke tests --- test/smoke-tests/circular-orbits.jl | 30 +++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 test/smoke-tests/circular-orbits.jl diff --git a/test/smoke-tests/circular-orbits.jl b/test/smoke-tests/circular-orbits.jl new file mode 100644 index 00000000..b17a1fe5 --- /dev/null +++ b/test/smoke-tests/circular-orbits.jl @@ -0,0 +1,30 @@ +# smoke test to make sure circular orbits work +using Test, Gradus + +@testset "pointfunctions" begin + + @testset "solve_equitorial_circular_orbit" begin + # only implemented for the BoyerLindquist metrics at the moment + # expected is sum of the circular orbit vϕ + # last updated: 22 Apr 2022 + r_range = 6.0:0.5:10.0 + for (m, expected) in [ + (BoyerLindquistAD(M = 1.0, a = 0.0), 0.5432533297869712), + (BoyerLindquistAD(M = 1.0, a = 1.0), 0.5016710246454921), + (BoyerLindquistAD(M = 1.0, a = -1.0), 0.5993458160081419), + (JohannsenAD(M = 1.0, a = 1.0, α22 = 1.0), 0.4980454719932759), + ] + vϕs = solve_equitorial_circular_orbit(m, r_range) + @test isapprox(sum(vϕs), expected, atol = 1e-6, rtol = 0.0) + end + end + + @testset "trace_equitorial_circular_orbit" begin + m = BoyerLindquistAD(M = 1.0, a = 0.0) + Gradus.isco(m) + r_range = 6.0:0.1:10.0 + simsols = @time trace_equitorial_circular_orbit(m, r_range) + # smoke test passed + @test true + end +end From 9fed8884a8e65d600dffe609f4a25bfb02b6ae86 Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Fri, 22 Apr 2022 19:30:27 +0100 Subject: [PATCH 8/8] renamed test-set --- test/smoke-tests/circular-orbits.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/smoke-tests/circular-orbits.jl b/test/smoke-tests/circular-orbits.jl index b17a1fe5..d3ac73de 100644 --- a/test/smoke-tests/circular-orbits.jl +++ b/test/smoke-tests/circular-orbits.jl @@ -1,7 +1,7 @@ # smoke test to make sure circular orbits work using Test, Gradus -@testset "pointfunctions" begin +@testset "circular-orbits" begin @testset "solve_equitorial_circular_orbit" begin # only implemented for the BoyerLindquist metrics at the moment