clean up

.gitignore

@@ -5,3 +5,4 @@
 Manifest.toml
 notes.md
 test/ad.jl
+test/mwe.jl

Project.toml

@@ -4,18 +4,23 @@ authors = ["Tom Kimpson"]
 version = "0.1.0"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extras]
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Test", "DoubleFloats"]

src/RelativisticDynamics.jl

@@ -6,7 +6,7 @@ module RelativisticDynamics
 # Imports
 import Parameters: @with_kw, @unpack
 import DifferentialEquations
-using Tullio,Combinatorics
+using Tullio,Combinatorics,LinearAlgebra
 
 # Exports
 export  orbit, PlotTrajectory, StackedPlot
@@ -16,7 +16,6 @@ include("system_parameters.jl")
 include("useful_functions.jl")
 include("metric.jl")
 include("universal_constants.jl")
-#include("model.jl")
 include("initial_conditions.jl")
 include("timestepping.jl")
 include("orbit.jl")

src/initial_conditions.jl

@@ -52,7 +52,6 @@ function initial_conditions(M)
 
     #Turn 4 velocity into 4 momentum
     pvector = m0*[tdot,rdot,θdot,ϕdot]
-
     pvector_covar = convert_to_covariant(g,pvector)
 
 

src/metric.jl

@@ -1,6 +1,5 @@
-#using Zygote
-#using PreallocationTools
-
+using SciMLSensitivity
+using Zygote
 
 """
     Δ = delta(r,a)
@@ -26,10 +25,12 @@ Construct the NxN matrix of the covariant metric.
 function covariant_metric(coords,a)
 
     g = zeros(typeof(a),4,4)
+
     t,r,θ,ϕ =  coords[1],coords[2],coords[3],coords[4]
     Σ = sigma(r,θ,a)
     Δ = delta(r,a)
 
+
     g[1,1] =   -(1.0 - 2.0*r / Σ)
     g[2,2] =   Σ / Δ
     g[3,3] =   Σ 

src/orbit.jl

@@ -8,6 +8,8 @@ function orbit(::Type{NF}=Float64;              # number format, use Float64 as
                kwargs...                        # all additional non-default parameters
                ) where {NF<:AbstractFloat}
 
+
+
     # Setup all system parameters, universal constants etc.
     P = SystemParameters(NF=NF;kwargs...) # Parameters
     bounds_checks(P)                      # Check all parameters are reasonable
@@ -16,9 +18,10 @@ function orbit(::Type{NF}=Float64;              # number format, use Float64 as
 
     #Initial conditions 
     initialization = initial_conditions(M)
-
+    
     #Evolve in time
     solution = timestepping(initialization, M)
+
     return solution, M
 
 end

src/plotting.jl

@@ -3,7 +3,6 @@
 
 using Plots 
 using LaTeXStrings
-#using Printf
 
 """
     PlotTrajectory(solution,model,dimensions=[1,2,3],savepath="")
@@ -15,8 +14,8 @@ function PlotTrajectory(solution,model,dimensions=[1,2,3],savepath="")
 
     @unpack a = model.parameters    #Get the BH spin parameter 
 
-    println("Plotting the solution generated with the following user defined parameters")
-    display(model.parameters)
+    println("Plotting the solution generated with the below user-defined parameters")
+    println(model.parameters)
     println("-------------------------------")
 
     #Interpolate to higher resolution for smooth plotting   

src/timestepping.jl

@@ -1,7 +1,5 @@
 
 using DifferentialEquations
-#using ComponentArrays
-#using SciMLSensitivity
 using Parameters: @unpack
 
 
@@ -11,20 +9,15 @@ The timestepping integration once all variables have been initialised
 """
 function timestepping(X::PrognosticVariables, M)
 
-
 @unpack a,m0, Tint = M.constants
 
-
 tspan = (zero(M.parameters.NF),M.constants.Tint) 
 u = vcat(X.xvector,X.pvector,X.svector)
 params = [a,m0]
 
 
 ode_prob = ODEProblem(MPD!,u,tspan,params)
 ode_solution = solve(ode_prob,DifferentialEquations.RK4())
-#ode_solution = solve(ode_prob,DifferentialEquations.Tsit5()) #, reltol=1e-8, abstol=1e-8)
-
-
 return ode_solution
 
 

src/universal_constants.jl

@@ -48,7 +48,7 @@ function Constants(P::SystemParameters)
     #Initial coordinates
     @unpack α = P
     r_initial = α
-    θ_initial = π/2.0
+    θ_initial = π/2.0 # starts in the plane by default
     ϕ_initial = 0.0
 
 
@@ -70,9 +70,6 @@ function Constants(P::SystemParameters)
     @unpack Norbits = P
     Tint = Norbits*2*π*α^(3/2)
 
-
-
-
 
    # This implies conversion to NF
     return Constants{P.NF}(light_c,μ,
@@ -83,7 +80,6 @@ function Constants(P::SystemParameters)
                            P.a)
 
 
-
 end
 
 

src/useful_functions.jl

@@ -1,7 +1,3 @@
-using LinearAlgebra
-#using ChainRulesCore
-
-
 """
     p_{μ} = convert_to_covariant(metric,p^{μ})
 Convert a vector from contravariant form to convariant form using the covariant metric 
@@ -23,8 +19,6 @@ function levi_civita_symbol()
 
     # Levi civita tensor
     levi = zeros(Integer,4,4,4,4) 
-    #ChainRulesCore.ignore_derivatives() do # This can be safely ignored by the differentiator - no dependence on the input parameters.
-
         for i in 1:4
             for j in 1:4
                 for k in 1:4
@@ -34,9 +28,7 @@ function levi_civita_symbol()
                     end
                 end 
             end
-        end 
-
-    #end  
+        end  
 
     return levi 
 end 

test/initial_conditions.jl

@@ -1,12 +1,3 @@
-
-using RelativisticDynamics
-using Test
-#using Zygote
-#using LinearAlgebra
-using Distributions
-using Tullio
-
-
 @testset "Initial conditions runs OK for default initialization" begin
 
     NF = Float64

test/metric.jl

@@ -1,9 +1,3 @@
-using RelativisticDynamics
-using Test
-using Enzyme
-using LinearAlgebra
-using Distributions
-using Tullio
 
 @testset "Trace of the metric" begin
 

test/number_formats.jl

@@ -1,12 +1,4 @@
 
-
-
-using RelativisticDynamics
-using Test
-import Parameters: @with_kw, @unpack
-using Tullio
-
-
 @testset "Check all number formats are consistent" begin
 
     for NF in [Float64,Float32]

test/orbit.jl

@@ -1,7 +1,3 @@
-using RelativisticDynamics
-using Test
-
-
 @testset "Basic run through" begin
     try
         solution,model = RelativisticDynamics.orbit()

test/runtests.jl

@@ -1,8 +1,9 @@
 using RelativisticDynamics
 using Test
-#using Zygote
-#using LinearAlgebra
+using LinearAlgebra
 using Distributions
+using Tullio
+using Enzyme
 
 # GENERAL
 include("metric.jl")

test/universal_constants.jl

@@ -1,17 +1,3 @@
-using RelativisticDynamics
-using Test
-#using Zygote
-#using LinearAlgebra
-using Distributions
-
-
-using Combinatorics
-
-
-
-
-
-
 @testset "MPD Carter constant for zero inclination" begin
 
     NF = Float64

test/useful_functions.jl

@@ -1,9 +1,6 @@
-
-
-
 @testset "Raise and lower vector indices" begin
 
-    for n in 1:1
+    for n in 1:5
         r = rand(Uniform(3.0,1e5))      # Radial coordinate. 3.0 as rough lower limit of an event horizon
         θ = rand(Uniform(0.0, 2.0*π))
         a = rand(Uniform(0.0, 1.0))
@@ -24,6 +21,27 @@
 end
 
 
+@testset "Levi civita is antisymmetric" begin
+
+    for n in 1:5
+        r = rand(Uniform(3.0,1e5))      # Radial coordinate. 3.0 as rough lower limit of an event horizon
+        θ = rand(Uniform(0.0, 2.0*π))
+        a = rand(Uniform(0.0, 1.0))
+
+
+        xvector=[0.0,r,θ,0.0] 
+        g  = RelativisticDynamics.covariant_metric(xvector,a)
+
+        levi = RelativisticDynamics.levi_civita_tensor(g)  #This is the fully contravariant Levi Civita tensor 
+        @test isapprox(permutedims(levi, (2,1,3,4)),-levi) #exchange of just 1 and 2
+
+    end 
+
+end
+
+
+
+