Merge pull request #42 from tomkimpson/langfzac_comments_1

langfzac comments

.github/workflows/CI.yml

@@ -18,7 +18,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.7'
+          - '1.9' #Latest
+          - '1.6' #LTS
         os:
           - ubuntu-latest
         arch:

Project.toml

@@ -4,23 +4,20 @@ 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"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 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"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "DoubleFloats"]
+test = ["Distributions", "Documenter", "Test"]

docs/src/how_to_run.md

@@ -9,19 +9,28 @@ solution,model = orbit()
 
 The `orbit()` function returns two objects. The first, `solution` holds the evolution of the position, momentum and spin vectors. The second, `model`, holds a copy of all the parameters and settings used to generate the solution, e.g. what was the BH spin?
 
-All default parameters can be found in `src/default_parameters.jl`. Passing a keyword argument to `orbit()` overrides the defaults e.g.
+All default parameters can be found in `src/system_parameters.jl`. Passing a keyword argument to `orbit()` overrides the defaults e.g.
 
-```julia
-run_speedy(e=0.6,a=0.99)
+```jldoctest
+using RelativisticDynamics
+solution,model = orbit(e=0.6,a=0.99);
+model.constants.a
+# output
+0.99
 ```
-would generate the solution for a system with an eccentricity = 0.6, around a BH with an extremal spin. 
+generates the solution for a system with an eccentricity = 0.6, around a BH with an extremal spin. 
 
 If provided, the number format has to be the first argument, all other arguments are keyword arguments. e.g. 
-```julia
-run_speedy(Float32,e=0.6,a=0.99)
+
+```jldoctest
+using RelativisticDynamics
+solution,model = orbit(Float32,e=0.6,a=0.99);
+model.parameters.NF
+# output
+Float32
 ```
 
-Please see `notebooks/demo.ipynb` for some worked examples using `RelativisticDynamics.jl`, including the application of autodiff methods.
 
+Please see `notebooks/demo.ipynb` for some worked examples using `RelativisticDynamics.jl`, including the application of autodiff methods.
 
 

docs/src/index.md

@@ -41,11 +41,16 @@ Direct installation from the git branch is also possible:
 ```
 
 
+## Troubleshooting and contributions 
+
+For any questions or problems, please feel free to [open an issue](https://github.com/tomkimpson/RelativisticDynamics.jl/issues)
+
+Contributions are always welcome - just [open a pull request](https://github.com/tomkimpson/RelativisticDynamics.jl/pulls).
+
 
 ## Acknowledgements 
 This work solving the MPD equations was originally motivated through interesting discussions with [Kinwah Wu](https://www.ucl.ac.uk/mssl/people/prof-kinwah-wu). The port to a modern, precision-flexible model in Julia was heavily inspired by [Milan Klöwer](https://github.com/milankl). Huge thanks to both.
 
 
-Contributions are always welcome - just [open a pull request](https://github.com/tomkimpson/RelativisticDynamics.jl/pulls)
 
 

docs/src/visualisation.md

@@ -5,19 +5,15 @@ It is often desireable, as a sanity check, to plot the orbital solution. Once th
 using Plots
 using RelativisticDynamics
 
-solution,model = orbit()  
-plot(solution,idxs=[1,2]) # to plot the time evolution of the first and second variables, in this case t and r 
+solution,model = orbit()
+plot(solution,model.constants.a) 
 ```
-
-Alternatively the user can use the inbuilt `PlotTrajectory()` function to plot the spatial path in usual Cartesian `x,y,z` coords, in either 3D or 2D. Here the indexes 1,2,3 refer to the `x`,`y`,`z` coordinates respectively. 
+By default the plotting recipe plots the $x-y$ orbital trajectory and up-samples the numerical solution by a factor of $10$. The user can define their own up-sampling factor, and pass any two of the integration variables as arguments, as well as the Cartesian coordinates ($x,y,z$). For example:
 
 ```julia
-PlotTrajectory(solution,model,[1,2],"example_media/2D_example.png")   # Plot a 2D example, save to example_media/2D_example.png
-PlotTrajectory(solution,model,[1,2,3],"example_media/3D_example.png") # Plot a 3D example, save to example_media/2D_example.png
-```
-
+plot(solution,model.constants.a,upsample=2, vars_to_plot = [:t,:r])     # Plot a timeseries of the r-coordinate, upsampled by a factor of 2
+plot(solution,model.constants.a,upsample=100, vars_to_plot = [:sθ,:sϕ]) # Plot the θ-ϕ components of the spin vector against each other, upsampled by a factor of 100
+plot(solution,model.constants.a,vars_to_plot = [:x,:z]) # Plot the x-z orbital trajectory
 
-One can also create a stacked plot of the trajectory in the $x-y$ and $x-z$ planes:
-```julia
-StackedPlot(solution,model,"../example_media/e08_stacked.pdf")
 ```
+

paper/paper.bib

@@ -168,3 +168,31 @@ @ARTICLE{Witzany2022
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+@article{Andrade2021, 
+         doi = {10.21105/joss.03703}, 
+         url = {https://doi.org/10.21105/joss.03703}, 
+         year = {2021}, 
+         publisher = {The Open Journal}, 
+         volume = {6}, 
+         number = {68}, 
+         pages = {3703}, 
+         author = {Tomas Andrade and Llibert Areste Salo and Josu C. Aurrekoetxea and Jamie Bamber and Katy Clough and Robin Croft and Eloy de Jong and Amelia Drew and Alejandro Duran and Pedro G. Ferreira and Pau Figueras and Hal Finkel and Tiago Fran\c{c}a and Bo-Xuan Ge and Chenxia Gu and Thomas Helfer and Juha Jäykkä and Cristian Joana and Markus Kunesch and Kacper Kornet and Eugene A. Lim and Francesco Muia and Zainab Nazari and Miren Radia and Justin Ripley and Paul Shellard and Ulrich Sperhake and Dina Traykova and Saran Tunyasuvunakool and Zipeng Wang and James Y. Widdicombe and Kaze Wong}, 
+         title = {GRChombo: An adaptable numerical relativity code for fundamental physics}, 
+         journal = {Journal of Open Source Software} }
+
+
+
+@article{PhysRevD.45.1840,
+  title = {Strong-field tests of relativistic gravity and binary pulsars},
+  author = {Damour, Thibault and Taylor, J. H.},
+  journal = {Phys. Rev. D},
+  volume = {45},
+  issue = {6},
+  pages = {1840--1868},
+  numpages = {0},
+  year = {1992},
+  month = {Mar},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevD.45.1840},
+  url = {https://link.aps.org/doi/10.1103/PhysRevD.45.1840}
+}

paper/paper.md

@@ -57,7 +57,7 @@ Together, equations \ref{eq:mpd1} - \ref{eq:mpd3} form the Mathisson-Papetrou-Di
 <!-- A Statement of need section that clearly illustrates the research purpose of the software and places it in the context of related work -->
 # Statement of need
 
-`RelativisticDynamics.jl` is an open-source Julia package for relativistic spin-orbital dynamics in the gravitational strong field for a Kerr spacetime. Existing codes for modelling the dynamics of spinning objects like pulsars in the strong-field regime are generally lacking, since such systems occupy an intermediate regime that is generally overlooked. At the "low" end of this regime there are post-Newtonian or geodesic descriptions which neglect the influence of the pulsar spin on the underlying spacetime metric ("spin-curvature" coupling). At the "high" end there is the full Numerical Relativity (NR) solutions which are primarily applicable to two BHs with a mass ratio $\mathcal{O}(1)$, and are computationally intractable for these MSP systems which are observed over a large number of orbital cycles.  
+`RelativisticDynamics.jl` is an open-source Julia package for relativistic spin-orbital dynamics in the gravitational strong field for a Kerr spacetime. Existing codes for modelling the dynamics of spinning objects like pulsars in the strong-field regime are generally lacking, since such systems occupy an intermediate regime that is generally overlooked. At the "low" end of this regime there are post-Newtonian or geodesic descriptions [e.g. @PhysRevD.45.1840] which neglect the influence of the pulsar spin on the underlying spacetime metric ("spin-curvature" coupling). At the "high" end there is the full Numerical Relativity (NR) solutions [e.g. @Andrade2021] which are primarily applicable to two BHs with a mass ratio $\mathcal{O}(1)$, and are computationally intractable for these MSP systems which are observed over a large number of orbital cycles.  
 
 `RelativisticDynamics.jl` aims to bridge this gap by providing a modern, fast code for accurate numerical evolution of spinning relativistic systems, via the MPD formalism. Julia is a modern language that solves the "two language problem", enabling fast dynamic typing and JIT compilation in conjunction with petaflop performance, comparable with numerical languages that are better known in the astrophysics community such as C or Fortran. As a modern language, it also provides a dedicated package manager and a large catalogue of _composable_ packages for scientific computing. This enables `RelativisticDynamics.jl` to easily leverage and interface with other scientific computing packages. The author and collaborators have used the general methods and mathematics described in this package for multiple research projects [e.g. @Kimpson2019; @Li2019; @Kimpson2020; @KimpsonAA] with a particular focus on the radio signals from spinning pulsar systems.  This package represents an attempt to create a documented, well-tested, open source resource for public use in this area, that can also be used as a computational playground for exploring techniques that could be applicable to more advanced numerical models. The package has been formulated in terms of ODE integration, rather than using e.g. action-angle variables [@Witzany2022], to allow for extension to general spacetime metrics and straightforward computatin of quantities relevant for pulsar observations e.g. spin axis orientation. 
 

src/RelativisticDynamics.jl

@@ -1,15 +1,13 @@
 module RelativisticDynamics
 
-
-
-
 # Imports
 import Parameters: @with_kw, @unpack
 import DifferentialEquations
 using Tullio,Combinatorics,LinearAlgebra
+using Zygote
 
 # Exports
-export  orbit, PlotTrajectory, StackedPlot
+export  orbit
 
 #Includes
 include("system_parameters.jl")

src/initial_conditions.jl

@@ -9,15 +9,13 @@ struct PrognosticVariables{NF<:AbstractFloat}
     svector         ::AbstractVector{NF}       
 end
 
-
 """
     initialization = initial_conditions(M)
 Setup the initial conditions for the MPD orbital dynamics"""
 function initial_conditions(M)
 
-
     @unpack NF = M.parameters
-   
+
     # 1. Four- position
     @unpack r_initial,θ_initial, ϕ_initial = M.constants
     xvector = [0.0,r_initial,θ_initial, ϕ_initial] # Starting coordinates
@@ -27,13 +25,11 @@ function initial_conditions(M)
     r,θ= xvector[2],xvector[3] 
     Δ = delta(r,a)
     Σ = sigma(r,θ,a)
-    g = covariant_metric(xvector,a)
-
+    g= covariant_metric(xvector,a) 
 
     # 2. Four - momentum 
     @unpack E,L,Q,m0 = M.constants 
 
-
     #These are 4 velocities from Schmidt 2002.
     #Initial Rdot is +ve as standard
     Pbar = E*(r^2+a^2) - a*L 
@@ -42,20 +38,15 @@ function initial_conditions(M)
     θbar = Q - ((1-E^2)*a^2 + L^2/sin(θ)^2)*cos(θ)^2
     ϕbar = a*Pbar/Δ -a*E + L/sin(θ)^2
 
-
     tdot = Tbar/Σ
     rdot = sqrt(Rbar)/Σ
     θdot = sqrt(θbar)/Σ
     ϕdot = ϕbar/Σ
 
-
-
     #Turn 4 velocity into 4 momentum
     pvector = m0*[tdot,rdot,θdot,ϕdot]
     pvector_covar = convert_to_covariant(g,pvector)
 
-
-
     # #3. Four - spin
     @unpack Sθ,Sϕ = M.parameters
     @unpack s0 = M.constants 

src/metric.jl

@@ -1,6 +1,3 @@
-using SciMLSensitivity
-using Zygote
-
 """
     Δ = delta(r,a)
 The well-known delta function of the Kerr metric
@@ -18,39 +15,72 @@ return r^2 + a^2 * cos(θ)^2
 end 
 
 
-"""
-    g=covariant_metric(coords,a)
-Construct the NxN matrix of the covariant metric.
-"""
+# """
+#     g=covariant_metric(coords,a)
+# Construct the NxN matrix of the covariant metric.
+# """
+# function covariant_metric(coords,a)
+
+#     #xs = zeros(typeof(a),4,4)
+#     #g = Zygote.Buffer(xs) #zeros(typeof(a),4,4)
+
+#     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] =   Σ 
+#     g[4,4] =   sin(θ)^2 * ((r^2 +a^2)^2 - Δ*a^2*sin(θ)^2) / Σ
+#     g[1,4] =   -2.0*a*r*sin(θ)^2/Σ 
+#     g[4,1] =   g[1,4] 
+
+#    # println("Normal metric")
+#    # display(g)
+
+#     #return copy(g)
+#     return g
+
+# end 
+
+
+
 function covariant_metric(coords,a)
 
-    g = zeros(typeof(a),4,4)
+    xs = zeros(typeof(a),4,4)
+    g = Zygote.bufferfrom(xs) 
 
     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] =   Σ 
     g[4,4] =   sin(θ)^2 * ((r^2 +a^2)^2 - Δ*a^2*sin(θ)^2) / Σ
     g[1,4] =   -2.0*a*r*sin(θ)^2/Σ 
     g[4,1] =   g[1,4] 
 
-    return g
+    return copy(g)
 
 end 
 
 
+
 """
     g=contravariant_metric(coords,a)
 Construct the NxN matrix of the contravariant metric.
 Metric components are defined via indvidual functions to allow for auto diff in unit tests
 """
 function contravariant_metric(coords,a)
 
-    metric_contra = zeros(typeof(a),4,4)
+
+
+    xs = zeros(typeof(a),4,4)
+    metric_contra = Zygote.bufferfrom(xs)
     t,r,θ,ϕ =  coords[1],coords[2],coords[3],coords[4]
     Σ = sigma(r,θ,a)
     Δ = delta(r,a)
@@ -64,10 +94,12 @@ function contravariant_metric(coords,a)
     metric_contra[4,4] = -(-(1.0 - 2.0*r / Σ))/denom
     metric_contra[1,4] = (-2.0*a*r*sin(θ)^2/Σ)/denom
     metric_contra[4,1] = metric_contra[1,4]
-    return metric_contra
+    return copy(metric_contra)
 end 
 
 
+
+
 """
     Γ = christoffel(coords,a)
 The christoffel symbols of the Kerr metric. See e.g. https://arxiv.org/pdf/0904.4184.pdf

src/orbit.jl

@@ -1,4 +1,3 @@
-
 """
     solution,model = orbit(NF,kwargs...)
     
@@ -7,21 +6,18 @@ Runs RelativisticDynamics.jl with number format `NF` and any additional paramete
 function orbit(::Type{NF}=Float64;              # number format, use Float64 as default
                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
     C = Constants(P)                      # Constants
     M = Model(P,C)                        # Pack all of the above into a single *Model struct 
 
-    #Initial conditions 
+    #Initial conditions
     initialization = initial_conditions(M)
 
     #Evolve in time
     solution = timestepping(initialization, M)
-    
+
     return solution, M
 
 end