From 483b701ccfdb5ec51f012f12007975a33ed39c7c Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Mon, 14 Oct 2024 17:00:11 +0100 Subject: [PATCH 1/6] fix: EMDA metric, correct expressions and inner radius Calculated the inner radius by hand, as there is currently no good cache support for this value, so having a root solver here would be a considerable performance overhead. Also re-implemented the metric and removed the fastmath macros to give better stability and remove ambiguous UTF-8 characters. --- src/metrics/dilaton-axion-ad.jl | 54 ++++++++++++++++----------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/src/metrics/dilaton-axion-ad.jl b/src/metrics/dilaton-axion-ad.jl index 64876a32..6bacbbbb 100644 --- a/src/metrics/dilaton-axion-ad.jl +++ b/src/metrics/dilaton-axion-ad.jl @@ -5,42 +5,41 @@ García et al. (1995). module __DilatonAxionAD using ..StaticArrays -@fastmath Σ(r, a, θ) = r^2 + a^2 * cos(θ)^2 -@fastmath Δ(r, a, rg) = r^2 - 2rg * r + a^2 -@fastmath A(r, a, Δ, θ) = (r^2 + a^2)^2 - a^2 * Δ * sin(θ)^2 +Σ(r, a, θ) = r^2 + a^2 * cos(θ)^2 +Δ(r, R, a) = r^2 + a^2 - R * r +# the a here is actually the dimensionless spin parameter a / R ? +Σhat(Σ, β, b, r, βb, a, θ, R) = Σ - (β^2 + 2b * r) + R^2 * βb * (βb - 2 * a * cos(θ)) +Δhat(Δ, β, b, r, βb, R) = Δ - (β^2 + 2b * r) - R * (R + 2b) * βb^2 -@fastmath W(βab, θ, βa) = 1 + (βab * (2 * cos(θ) - βab) + βa^2) * csc(θ)^2 -@fastmath Σ̂(Σ, β, b, r, βb, a, θ, rg) = - Σ - (β^2 + 2b * r) + rg^2 * βb * (βb - 2 * (a / rg) * cos(θ)) -@fastmath Δ̂(Δ, β, b, r, βb, rg) = Δ - (β^2 + 2b * r) - rg * (rg + 2b) * βb^2 -@fastmath Â(δ, a, Δ̂, W, θ) = δ^2 - a^2 * Δ̂ * W^2 * sin(θ)^2 -@fastmath δ(r, b, a) = r^2 - 2b * r + a^2 +δ(r, b, a) = r^2 - 2b * r + a^2 + +W(θ, βab, βa) = 1 + (βab * (2 * cos(θ) - βab) + βa^2) * csc(θ)^2 +A(a, θ, δ, W, Δhat) = δ^2 - Δhat * (W * a * sin(θ))^2 @fastmath function metric_components(M, a, β, b, rθ) (r, θ) = rθ - rg = M - - Δ₀ = Δ(r, a, rg) - Σ₀ = Σ(r, a, θ) + R = M - βa = β / a + # check these for the normalisation term R βb = β / b - βab = rg * β / (a * b) + βa = β / a + βab = β / (a * b) - Δ̂₀ = Δ̂(Δ₀, β, b, r, βb, rg) - Σ̂₀ = Σ̂(Σ₀, β, b, r, βb, a, θ, rg) + Σ₀ = Σ(r, a, θ) + Δ₀ = Δ(r, R, a) + Δhat₀ = Δhat(Δ₀, β, b, r, βb, R) + Σhat₀ = Σhat(Σ₀, β, b, r, βb, a, θ, R) δ₀ = δ(r, b, a) - W₀ = W(βab, θ, βa) - - Â₀ = Â(δ₀, a, Δ̂₀, W₀, θ) + W₀ = W(θ, βab, βa) + A₀ = A(a, θ, δ₀, W₀, Δhat₀) - tt = -(Δ̂₀ - a^2 * sin(θ)^2) / Σ̂₀ - rr = Σ̂₀ / Δ̂₀ - θθ = Σ̂₀ - ϕϕ = Â₀ * sin(θ)^2 / Σ̂₀ + tt = -(Δhat₀ - a^2 * sin(θ)^2) / Σhat₀ + rr = Σhat₀ / Δhat₀ + θθ = Σhat₀ + ϕϕ = A₀ * sin(θ)^2 / Σhat₀ - tϕ = -a * (δ₀ - Δ̂₀ * W₀) * sin(θ)^2 / Σ̂₀ + tϕ = -a * (δ₀ - Δhat₀ * W₀) * sin(θ)^2 / Σhat₀ @SVector [tt, rr, θθ, ϕϕ, tϕ] end @@ -60,7 +59,7 @@ Einstein-Maxwell-Dilaton-Axion metric. "Singularity mass." M = 1.0 "Singularity spin." - a = 0.0 + a = 0.5 "Dilaton coupling strength." β = 0.0 "Axion coupling strength." @@ -69,6 +68,7 @@ end metric_components(m::DilatonAxion{T}, rθ) where {T} = __DilatonAxionAD.metric_components(m.M, m.a, m.β, m.b, rθ) -inner_radius(m::DilatonAxion{T}) where {T} = m.M + √(m.M^2 - m.a^2) +inner_radius(m::DilatonAxion{T}) where {T} = + m.M + m.b + √((m.M + m.b)^2 - m.a^2 + m.β^2 - (m.M - 2m.b) * m.M * (m.β / m.b)^2) export DilatonAxion From ab7b7bbf96c7cd60f9447149f274ef8111f358c8 Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Mon, 14 Oct 2024 17:01:31 +0100 Subject: [PATCH 2/6] feat: velocity interp now stores polar coordinate We now also keep the theta (polar) velocity array and return it when interpolating the plunging region velocities. This is because, e.g. the EMDA metric in free-fall has a non-zero polar component, which needs to be taken account of. --- src/orbits/orbit-solving.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/orbits/orbit-solving.jl b/src/orbits/orbit-solving.jl index 0c91e9fd..00a5a560 100644 --- a/src/orbits/orbit-solving.jl +++ b/src/orbits/orbit-solving.jl @@ -100,6 +100,7 @@ struct PlungingInterpolation{M,_interp_type} m::M t::_interp_type r::_interp_type + θ::_interp_type ϕ::_interp_type function PlungingInterpolation(m::M, sol) where {M<:AbstractMetric{T}} where {T} @@ -110,6 +111,7 @@ struct PlungingInterpolation{M,_interp_type} vt = sol[5, :][I] vr = sol[6, :][I] + vθ = sol[7, :][I] vϕ = sol[8, :][I] r_interp = _make_interpolation(r, vt) @@ -117,6 +119,7 @@ struct PlungingInterpolation{M,_interp_type} m, r_interp, _make_interpolation(r, vr), + _make_interpolation(r, vθ), _make_interpolation(r, vϕ), ) end @@ -130,8 +133,9 @@ function (pintrp::PlungingInterpolation)(r) r_bounded = _enforce_interpolation_bounds(r, pintrp) vt = pintrp.t(r_bounded) vr = pintrp.r(r_bounded) + vθ = pintrp.r(r_bounded) vϕ = pintrp.ϕ(r_bounded) - SVector(vt, vr, 0, vϕ) + SVector(vt, vr, vθ, vϕ) end function interpolate_plunging_velocities( From 360a1712317cd622a30822d20b84f1fb7fc6726f Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Mon, 14 Oct 2024 17:02:58 +0100 Subject: [PATCH 3/6] feat: allow offset as kwarg in plunging interp --- src/orbits/orbit-solving.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/orbits/orbit-solving.jl b/src/orbits/orbit-solving.jl index 00a5a560..fad006e1 100644 --- a/src/orbits/orbit-solving.jl +++ b/src/orbits/orbit-solving.jl @@ -143,12 +143,12 @@ function interpolate_plunging_velocities( max_time = 50_000, contra_rotating = false, reltol = 1e-9, + δr = (reltol * 10), kwargs..., ) where {T} isco = Gradus.isco(m) # rule of thumb to achieve desired error - δr = (reltol * 10) u = @SVector([0.0, isco - δr, deg2rad(90), 0.0]) v = Gradus.CircularOrbits.plunging_fourvelocity( m, From 55e88b7a372cac3e793adaa98d87602a400480fe Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Mon, 14 Oct 2024 17:17:52 +0100 Subject: [PATCH 4/6] test: update the ISCO test for EMDA --- test/smoke-tests/special-radii.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/smoke-tests/special-radii.jl b/test/smoke-tests/special-radii.jl index 0d87158e..3d89ce83 100644 --- a/test/smoke-tests/special-radii.jl +++ b/test/smoke-tests/special-radii.jl @@ -12,7 +12,7 @@ using StaticArrays KerrMetric(1.0, -0.998), JohannsenMetric(M = 1.0, a = 0.998, α13 = 1.0), JohannsenMetric(M = 1.0, a = 0.998, α22 = 1.0), - DilatonAxion(M = 1.0, a = 0.998, β = 0.2, b = 1.0), + DilatonAxion(M = 1.0, a = 0.6, β = -0.5, b = 0.1), ) @testset "iscos" begin @@ -26,7 +26,7 @@ using StaticArrays 8.99437445480357, 2.8482863127671534, 1.1306596884484472, - 6.880202293032178, + 32.95283734688169, ], ) isco_r = Gradus.isco(m) From ba0e64e28951f2a42923ac80f0c40ee997b9987d Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Mon, 14 Oct 2024 17:25:45 +0100 Subject: [PATCH 5/6] test: update regression values For the EMDA, now testing against a value from a figure in Younsi et al, 2023. For the precisions, this I think is due to a change in Julia version. --- test/integration/test-precision.jl | 10 +++++----- test/transfer-functions/test-thick-disc.jl | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/test/integration/test-precision.jl b/test/integration/test-precision.jl index 8a5212b0..bc8bc4c4 100644 --- a/test/integration/test-precision.jl +++ b/test/integration/test-precision.jl @@ -6,12 +6,12 @@ m = KerrMetric(M = 1.0, a = 1.0) u = SVector(0.0, 1000.0, π / 2, 0.0) # target position -target = SVector(10.0, 0.001, 0.0) +target = SVector(10.0, 0.005, 0.0) α, β, accuracy = Gradus.impact_parameters_for_target(target, m, u) -@test α ≈ -0.0017523796800187875 rtol = 1e-3 +@test α ≈ -0.004013630261097743 rtol = 1e-3 @test β ≈ 10.969606493445841 rtol = 1e-3 -@test accuracy ≈ 0.009091572709970781 rtol = 1e-3 +@test accuracy ≈ 0.004587323209289997 rtol = 1e-3 # new target target = SVector(10.0, deg2rad(40), -π / 4) @@ -19,11 +19,11 @@ target = SVector(10.0, deg2rad(40), -π / 4) α, β, accuracy = Gradus.impact_parameters_for_target(target, m, u) @test α ≈ 4.848373364532467 rtol = 1e-3 @test β ≈ 8.012499652493656 rtol = 1e-3 -@test accuracy ≈ 0.0049975644770909105 rtol = 1e-3 +@test accuracy ≈ 0.0030663177403400794 rtol = 1e-3 # test geodesic point interface α, β, gp, accuracy = Gradus.optimize_for_target(target, m, u) @test gp.x[1] ≈ 1005.2700874611182 rtol = 1e-3 @test α ≈ 4.848373364532467 rtol = 1e-3 @test β ≈ 8.012499652493656 rtol = 1e-3 -@test accuracy ≈ 0.0049975644770909105 rtol = 1e-3 +@test accuracy ≈ 0.0030663177403400794 rtol = 1e-3 diff --git a/test/transfer-functions/test-thick-disc.jl b/test/transfer-functions/test-thick-disc.jl index b8637761..759e72f0 100644 --- a/test/transfer-functions/test-thick-disc.jl +++ b/test/transfer-functions/test-thick-disc.jl @@ -8,7 +8,7 @@ d = ShakuraSunyaev(m) tf = cunningham_transfer_function(m, x, d, 3.0; β₀ = 1.0) total = sum(filter(!isnan, tf.f)) -@test total ≈ 12.253422180875667 atol = 1e-4 +@test total ≈ 12.255956053708246 atol = 1e-4 m = KerrMetric(1.0, 0.2) x = SVector(0.0, 10_000, deg2rad(20), 0.0) @@ -16,7 +16,7 @@ d = ShakuraSunyaev(m; eddington_ratio = 0.2) tf = cunningham_transfer_function(m, x, d, 5.469668466100368; β₀ = 1.0) total = sum(filter(!isnan, tf.f)) -@test total ≈ 20.83469 atol = 1e-4 +@test total ≈ 20.83782174681574 atol = 1e-4 # the transfer function here is pretty horrible as it's almost impossible to actually # see; this is a test to make sure it doesn't error From 943709cb4a548ebdfac2033fa5f6b654502fa3e6 Mon Sep 17 00:00:00 2001 From: Fergus Baker Date: Mon, 14 Oct 2024 17:27:39 +0100 Subject: [PATCH 6/6] ci: bump julia version to 1.11 --- .github/workflows/auto-register.yml | 2 +- .github/workflows/docs.yml | 2 +- .github/workflows/test.yml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/auto-register.yml b/.github/workflows/auto-register.yml index 65f19c44..35d9e334 100644 --- a/.github/workflows/auto-register.yml +++ b/.github/workflows/auto-register.yml @@ -21,7 +21,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ['1.10'] + julia-version: ['1.11'] os: [ubuntu-latest] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index cedfbc42..7bc2b229 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -17,7 +17,7 @@ jobs: - uses: actions/checkout@v3 - uses: julia-actions/setup-julia@v1 with: - version: '1.10' + version: '1.11' - run: | julia --project=docs -e ' using Pkg diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index ead51ada..d6166b00 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -26,7 +26,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ['1.10'] + julia-version: ['1.11'] os: [ubuntu-latest] steps: