Refactor GRHD project

Major rework of GRHD project to include new broadcasting interface. The goal is to get a TOV evolution with the AV method.

PR contains various refactors:

• initialdata/: We promote this to a separate env in order to add CairoMakie for plotting.
• refactored and improve initialdata/SphericalAccretion.jl solver so that it agrees with various references: solver from bugner's thesis (extract from bamps), plots from two papers from Papadoulos.
• refactored initialdata/TOV.jl to not depend on dg1d.jl anymore, instead use OrindaryDiffEq and CairoMakie for plots
• copy of cons2prim from SRHD. these two versions should be merged eventually, but will do in a separate PR
• implement spherical1d formulation: this is a rewrite of the evolution GRHD equations to work with non-zero radial shift vector. this is needed to evolve the spherical accretion test in cowling approximation where there should be no dynamics.
• implement rescaled_spherical1d formulation: this is a rewrite of the previous formulation based on the generalized Valencia formulation (cf. Montero, Baumgarte paper from 2017); the auxiliary metric is chosen such that coordinate singularities at r=0 in the source terms are cancelled; doesn't work atm
• allow the spherical1d to include r=0 in the domain by using symmetry conditions to evaluate source terms at r=0 (workaround for rescaled_spherical1d not working right now)
• new ref tests: bondi accretion and TOV star
• implement a differentiate! method on Mesh in order to compute derivatives by hand, if needed

initialdata/SphericalAccretion.jl

+module SphericalAccretion
+
+
+
 using CairoMakie
+using CairoMakie.FileIO
+using DelimitedFiles
 using Printf
+using HDF5
+using Roots
+using UnPack
 
 
 """
@@ -7,53 +16,51 @@ Solution taken from Bugner [2017], section 4.2.
 
 Kerr-Schild coordinates allow use to go beyond horizon towards the singularity.
 However, it becomes more difficult for the Newton Rapshon solver to determine a solution
-once beyond the horizon. To this end, try increasing `N_rs` whenever you lower `rmin` and 
+once beyond the horizon. To this end, try increasing `N_rs` whenever you lower `rmin` and
 the solver aborts inside the horizon.
 Alternatively, one could also add an adaptive stepper.
 
 Note: We implicitly assume θ = π/2, because then sin(θ) = 1, although,
 the sines should cancel out everywhere.
 """
-function solve()
-
-  # parameters
-  M = 1.0
-  N_rs = 100000
-  rmin, rmax = 0.1*M, 200*M
-  rc, ρc, Γ = 200*M, 7e-4, 4/3
+function solve(; rmin, rmax, rc, ρc, Γ, M=1.0, N_rs=1000,
+                 urc_sign=1)
 
-  urc = sqrt(M / (2 * rc))
+  urc = urc_sign*sqrt(M /(2 * rc))
   # cannot find solutions within these regions, cf. Michel 1972
   @show urc
   @assert rc >= 6*M
-  @assert urc <= 1 / 3
+  @assert urc^2 <= 1 / 3
 
-  Vc2 = urc^2 / (1 - 3 * urc^2)
-  K = Vc2 * (Γ - 1) / ( (Γ - 1 - Vc2) * Γ * ρc^(Γ - 1) )
+  Vc = urc^2 / (1 - 3 * urc^2)
+  K = Vc * (Γ - 1) / ( (Γ - 1 - Vc) * Γ * ρc^(Γ - 1) )
+  @show K
 
   # here we set θ = π/2 so that sin(θ) = 1
-  C1 = rc^2 * ρc * urc
   hc = 1 + K * Γ * ρc^(Γ - 1) / (Γ - 1)
   u0c = sqrt(1 - 2 * M / rc + urc^2)
+  C1 = rc^2 * ρc * urc
+  C2 = rc^2 * ρc * hc * u0c * urc
   C3 = hc^2 * u0c^2
 
   function F(ur, r)
     ρ = C1 / (r^2 * ur)
     h = 1 + K * Γ * ρ^(Γ - 1) / (Γ - 1)
     u0 = sqrt(1 - 2 * M / r + ur^2)
-    u0 = sqrt(1 - 2 * M / r + ur^2)
     return h^2 * u0^2 - C3
   end
 
-  function dF(ur, r)
-    ρ = C1 / (r^2 * ur)
-    h = 1 + K * Γ * ρ^(Γ - 1) / (Γ - 1)
-    u0 = sqrt(1 - 2 * M / r + ur^2)
-    dρdur = - C1 / (r^2 * ur^2)
-    dhdρ = K * Γ * ρ^(Γ - 2)
-    du0dur = ur / u0
-    return 2 * (u0^2 * h * dhdρ * dρdur + h^2 * u0 * du0dur)
-  end
+  # manually compute derviative and
+  # use Roots.find_zero((F_at_r,dF_at_r), ur_guess, Roots.Newton()) below
+  # function dF(ur, r)
+  #   ρ = C1 / (r^2 * ur)
+  #   h = 1 + K * Γ * ρ^(Γ - 1) / (Γ - 1)
+  #   u0 = sqrt(1 - 2 * M / r + ur^2)
+  #   dρdur = - C1 / (r^2 * ur^2)
+  #   dhdρ = K * Γ * ρ^(Γ - 2)
+  #   du0dur = ur / u0
+  #   return 2 * (u0^2 * h * dhdρ * dρdur + h^2 * u0 * du0dur)
+  # end
 
   # radial domain over which we solve
   rs = collect(LinRange(rmin, rmax, N_rs))
@@ -65,52 +72,275 @@ function solve()
     # - use the previous result for the remaining radi
     # The following converges to the correct solution:
     if r > 2*M
-      ur_guess = sqrt(1 - 2 * M / r) # asymptotic value for ur
+      ur_guess = urc_sign*sqrt(1 - 2 * M / r) # = sqrt(1+gtt), =^= asymptotic value for ur when u0 -> 1
     else
       ur_guess = urs[idx+1] # use previous iterate where asymptotic value is not defined
     end
     F_at_r(ur) = F(ur, r)
-    dF_at_r(ur) = dF(ur, r)
+    # dF_at_r(ur) = dF(ur, r)
     ur = try
-      Roots.find_zero((F_at_r,dF_at_r), ur_guess, Roots.Newton())
+      # Roots.find_zero((F_at_r,dF_at_r), ur_guess, Roots.Newton())
+      Roots.find_zero(F_at_r, ur_guess)
     catch err
       println()
-      error("Could not determine solution for 'ur' at radius r = $r. " *
-            "Remaining error: $err")
+      @info("Could not determine solution for 'ur' at radius r = $r. ")
+      rethrow(err)
     end
     # @show (ur, r)
-    if ur < 0.0
-      @warn "Aborting due to ur becoming negative at r = $r." 
-      @info "$(N_rs - idx) points skipped ..."
-      break
-    end
+    # if ur < 0.0
+    #   @warn "Aborting due to ur becoming negative at r = $r."
+    #   @info "$(N_rs - idx) points skipped ..."
+    #   break
+    # end
     urs[idx] = ur
   end
 
   ρs = @. C1 / (rs^2 * urs)
   ps = @. K * ρs^Γ
   ϵs = @. ps / ((Γ - 1) * ρs)
-  βr = @. 2 * M / (rs + 2 * M)
+  βr = @. 2 * M / (rs + 2 * M) # = β^r
   α = @. sqrt( rs / (rs + 2 * M) )
-
-  gtt_uu = @. - (2 * M + rs) / rs
-  gtr_uu = @. 2 * M / rs
-  grr_uu = @. (rs - 2 * M) / rs
+  γrr   = @. 1 + 2*M/rs
+  γrr_uu = @. 1/γrr
+  gtt_uu = @. -1/α^2
+  gtr_uu = @. βr/α^2
+  grr_uu = @. γrr_uu - βr*βr/α^2
+  gtt   = @. -α^2 + γrr*βr*βr
+  gtr   = @. γrr * βr
+  grr   = @. γrr
   u0s_d = @. sqrt(1 - 2 * M / rs + urs^2)
   u0s_u = @. gtt_uu * u0s_d + gtr_uu / grr_uu * (urs - gtr_uu * u0s_d)
+  # u0s_u_v2 = @. (gtt_uu * u0s_d + gtr_uu*grr*urs)/(1-gtr*gtr_uu)
+  # @show(maximum(abs, u0s_u .- u0s_u_v2))
   vrs = @. urs / (α * u0s_u) + βr / α
+  vrs_d = @. γrr * vrs
+  hs = @. 1 + K * Γ * ρs^(Γ - 1) / (Γ - 1)
+  res_C1 = @. rs^2 * ρs * urs - C1
+  res_C2 = @. rs^2 * ρs * hs * u0s_d * urs - C2
+  res_C3 = @. hs^2 * u0s_d^2 - C3
+  vs = @. sqrt(γrr * vrs^2)
+
+  return (; rs, urs, vrs, ρs, ps, ϵs, res_C1, res_C2, res_C3, vs, vrs_d,
+            rmin, rmax, rc, Vc, ρc, Γ, M, N_rs, K)
+
+end
+
+
+function solve_ref()
+  # reference data from bugner's solver SphericalAccretion.c (extracted from bamps)
+  M = 1.0
+  rmin, rmax = 1.8*M, 200*M
+  rc, ρc, Γ = 200*M, 7e-4, 4/3
+  solve(; rmin, rmax, rc, ρc, Γ, M)
+end
+
+function solve_papadopoulus_1998()
+  # arXiv:gr-qc/9803087v1
+  # works with Eddington-Finkelstein coords (same form as Bugner; but Bugner calls them Kerr-Schild?)
+  M = 1.0
+  rmin, rmax = 0.5*M, 50*M
+  rc, ρc, Γ = 400*M, 1e-2, 4/3
+  solve(; rmin, rmax, rc, ρc, Γ, M)
+end
+
+function solve_papadopoulus_1999()
+  # arXiv:gr-qc/9902018v1
+  # works with ingoing Eddington-Finkelstein coords, but the equations appear to be the same
+  # compare with the 1998 paper, in these coordinates G(r) = 0
+  M = 1.0
+  rmin, rmax = 1.5*M, 30*M
+  rc, ρc, Γ = 200*M, 1e-2, 5/3
+  solve(; rmin, rmax, rc, ρc, Γ, M)
+end
+
+
+function save(fname, data)
+  # save
+  filename = joinpath(@__DIR__, fname)
+  h5open(filename, "w") do file
+    @info "Saving results to '$filename' ..."
+    g = create_group(file, "bondi_accretion")
+    HDF5.attributes(g)["rc"]   = data.rc
+    HDF5.attributes(g)["ρc"]   = data.ρc
+    HDF5.attributes(g)["Vc"]   = data.Vc
+    HDF5.attributes(g)["Γ"]    = data.Γ
+    HDF5.attributes(g)["M"]    = data.M
+    HDF5.attributes(g)["K"]    = data.K
+    HDF5.attributes(g)["rmin"] = data.rmin
+    HDF5.attributes(g)["rmax"] = data.rmax
+    g["r"]  = data.rs
+    g["ur"] = data.urs
+    g["vr"] = data.vrs
+    g["ρ"]  = data.ρs
+    g["p"]  = data.ps
+    g["ϵ"]  = data.ϵs
+  end
+  return
+end
+
+
+function plot(data)
+
+  @unpack M, ρs, ps, ϵs, vrs, rs, res_C1, res_C2, res_C3, vrs_d = data
 
   fig = Figure()
-  ax_ρ   = Axis(fig[1,1], xlabel="r", ylabel="ρ")
-  ax_p   = Axis(fig[1,2], xlabel="r", ylabel="p")
-  ax_ϵ   = Axis(fig[2,1], xlabel="r", ylabel="ϵ")
-  ax_v_r = Axis(fig[2,2], xlabel="r", ylabel="v_r")
+  ax_ρ      = Axis(fig[1,1], xlabel=L"r", ylabel=L"\rho")
+  ax_p      = Axis(fig[1,2], xlabel=L"r", ylabel=L"p")
+  ax_ϵ      = Axis(fig[2,1], xlabel=L"r", ylabel=L"\epsilon")
+  ax_v_r    = Axis(fig[2,2], xlabel=L"r", ylabel=L"v^r")
+  # ax_res_C1 = Axis(fig[3,1], xlabel=L"r", ylabel=L"R_{C_1}")
+  # ax_res_C2 = Axis(fig[3,2], xlabel=L"r", ylabel=L"R_{C_2}")
+  # ax_res_C3 = Axis(fig[4,1], xlabel=L"r", ylabel=L"R_{C_3}")
+  # ax_v_r_d  = Axis(fig[4,2], xlabel=L"r", ylabel=L"v_r")
+
+  rhor = 2*M # radius of horizon
 
   lines!(ax_ρ, rs, ρs)
   lines!(ax_p, rs, ps)
   lines!(ax_ϵ, rs, ϵs)
   lines!(ax_v_r, rs, vrs)
+  # lines!(ax_res_C1, rs, res_C1)
+  # lines!(ax_res_C2, rs, res_C2)
+  # lines!(ax_res_C3, rs, res_C3)
+  # lines!(ax_v_r_d, rs, vrs_d)
+  vlines!(ax_ρ, rhor, linestyle=:dash)
+  vlines!(ax_p, rhor, linestyle=:dash)
+  vlines!(ax_ϵ, rhor, linestyle=:dash)
+  vlines!(ax_v_r, rhor, linestyle=:dash)
+  # vlines!(ax_v_r_d, rhor, linestyle=:dash)
+
+  # xlims!(ax_ρ, (1.8,10.0))
+  # xlims!(ax_p, (1.8,10.0))
+  # xlims!(ax_ϵ, (1.8,10.0))
+  # xlims!(ax_v_r, (1.8,10.0))
 
   # use save("filename.pdf", fig) to save to pdf
   return fig
 end
+
+
+function plot_compare_ref()
+
+  mydata = solve_ref()
+  csv = readdlm(joinpath(@__DIR__, "bondi-acc-ref.txt"))
+  refdata = (; rs=csv[:,1],
+               ρs=csv[:,2],
+               ϵs=csv[:,3],
+               ps=csv[:,4],
+               vrs=csv[:,5],
+               urs=csv[:,6],
+               ut2=csv[:,7])
+
+  fig = Figure()
+  ax_ρ      = Axis(fig[1,1], xlabel=L"r", ylabel=L"\rho")
+  ax_p      = Axis(fig[1,2], xlabel=L"r", ylabel=L"p")
+  ax_ϵ      = Axis(fig[2,1], xlabel=L"r", ylabel=L"\epsilon")
+  ax_vr     = Axis(fig[2,2], xlabel=L"r", ylabel=L"v^r")
+  ax_ur     = Axis(fig[3,1], xlabel=L"r", ylabel=L"u^r")
+
+  lines!(ax_ρ,   mydata.rs,  mydata.ρs)
+  lines!(ax_p,   mydata.rs,  mydata.ps)
+  lines!(ax_ϵ,   mydata.rs,  mydata.ϵs)
+  lines!(ax_vr,  mydata.rs,  mydata.vrs)
+  lines!(ax_ur,  mydata.rs,  mydata.urs)
+  lines!(ax_ρ,   refdata.rs, refdata.ρs, linestyle=:dash)
+  lines!(ax_p,   refdata.rs, refdata.ps, linestyle=:dash)
+  lines!(ax_ϵ,   refdata.rs, refdata.ϵs, linestyle=:dash)
+  lines!(ax_vr,  refdata.rs, refdata.vrs, linestyle=:dash)
+  lines!(ax_ur,  refdata.rs, refdata.urs, linestyle=:dash)
+
+  rhor = 2*mydata.M
+  vlines!(ax_ρ, rhor, linestyle=:dash, color=:grey)
+  vlines!(ax_p, rhor, linestyle=:dash, color=:grey)
+  vlines!(ax_ϵ, rhor, linestyle=:dash, color=:grey)
+  vlines!(ax_vr, rhor, linestyle=:dash, color=:grey)
+  vlines!(ax_ur, rhor, linestyle=:dash, color=:grey)
+
+  return fig
+
+end
+
+
+function plot_compare_papadopoulus_1998()
+
+  mydata = solve_papadopoulus_1998()
+  img_p = load("bondi-acc_papadopoulos1998_p.png")
+  img_ρ = load("bondi-acc_papadopoulos1998_rho.png")
+  # v = sqrt(γ_rr) v^r; is called total velocity in the paper (hidden in text)
+  # the label in the plot is just 'Velocity'
+  img_v = load("bondi-acc_papadopoulos1998_v.png")
+
+  xmin, xmax = 0.0, 15.0
+  ymin_ρ, ymax_ρ = 0, 120.0
+  ymin_p, ymax_p = 0, 2.5
+  ymin_vr, ymax_vr = -0.3, -0.2
+
+  fig = Figure()
+  ax_p  = Axis(fig[1,1], xlabel=L"r", ylabel=L"p")
+  ax_ρ  = Axis(fig[2,1], xlabel=L"r", ylabel=L"\rho")
+  ax_vr = Axis(fig[3,1], xlabel=L"r", ylabel=L"v^r")
+
+  image!(ax_ρ, xmin..xmax, ymin_ρ..ymax_ρ, rotr90(img_ρ))
+  image!(ax_p, xmin..xmax, ymin_p..ymax_p, rotr90(img_p))
+  image!(ax_vr, xmin..xmax, ymin_vr..ymax_vr, rotr90(img_v))
+
+  lines!(ax_ρ,  mydata.rs, mydata.ρs, linestyle=:dash, color=:orange)
+  lines!(ax_p,  mydata.rs, mydata.ps, linestyle=:dash, color=:orange)
+  # don't understand why they plot with a minus here
+  # there seem to be ingoing and outgoing solutions (urc_sign parameter),
+  # but the value for v never fits
+  lines!(ax_vr, mydata.rs, -mydata.vs, linestyle=:dash, color=:orange)
+
+  xlims!(ax_ρ,   (xmin,xmax))
+  xlims!(ax_p,   (xmin,xmax))
+  xlims!(ax_vr,  (xmin,xmax))
+  ylims!(ax_ρ,   (ymin_ρ,ymax_ρ))
+  ylims!(ax_p,   (ymin_p,ymax_p))
+  ylims!(ax_vr,  (ymin_vr,ymax_vr))
+
+  return fig
+end
+
+
+function plot_compare_papadopoulus_1999()
+
+  mydata = solve_papadopoulus_1999()
+  img_ϵ = load("bondi-acc_papadopoulos1999_eps.png")
+  img_ρ = load("bondi-acc_papadopoulos1999_rho.png")
+  img_vr = load("bondi-acc_papadopoulos1999_vr.png")
+
+  xmin, xmax = 1.0, 10.0
+  ymin_ρ, ymax_ρ = 0, 10.0
+  ymin_ϵ, ymax_ϵ = 0, 0.25
+  ymin_vr, ymax_vr = -1.0, 0.0
+
+  fig = Figure()
+  ax_ϵ  = Axis(fig[1,1], xlabel=L"r", ylabel=L"\epsilon")
+  ax_ρ  = Axis(fig[2,1], xlabel=L"r", ylabel=L"\rho")
+  ax_vr = Axis(fig[3,1], xlabel=L"r", ylabel=L"u^r")
+
+  image!(ax_ρ, xmin..xmax, ymin_ρ..ymax_ρ, rotr90(img_ρ))
+  image!(ax_ϵ, xmin..xmax, ymin_ϵ..ymax_ϵ, rotr90(img_ϵ))
+  image!(ax_vr, xmin..xmax, ymin_vr..ymax_vr, rotr90(img_vr))
+
+  lines!(ax_ρ,  mydata.rs, mydata.ρs, linestyle=:dash, color=:orange)
+  lines!(ax_ϵ,  mydata.rs, mydata.ϵs, linestyle=:dash, color=:orange)
+  # don't understand why they plot with a minus here
+  # there seem to be ingoing and outgoing solutions (urc_sign parameter),
+  # but the value for ur never fits
+  lines!(ax_vr, mydata.rs, -mydata.urs, linestyle=:dash, color=:orange)
+  display(extrema(mydata.urs))
+
+  xlims!(ax_ρ,   (xmin,xmax))
+  xlims!(ax_ϵ,   (xmin,xmax))
+  xlims!(ax_vr,  (xmin,xmax))
+  ylims!(ax_ρ,   (ymin_ρ,ymax_ρ))
+  ylims!(ax_ϵ,   (ymin_ϵ,ymax_ϵ))
+  ylims!(ax_vr,  (ymin_vr,ymax_vr))
+
+  return fig
+end
+
+
+
+end # module