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

src/GRHD/cons2prim.jl

@@ -110,14 +110,10 @@ end
 
 
 @inline function cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr)
-  # @accepts D, Sr, τ, r
-  # @accepts γrr
 
   rcoord = r
 
   @unpack atm_density, atm_threshold, eos = equation
-  # atm_density = equation.atmosphere.ρ
-  # atm_threshold = equation.atmosphere.threshold
   ρmin = atm_density * atm_threshold
   @unpack K, Gamma = eos
   pmin = K * ρmin^Gamma
@@ -126,9 +122,7 @@ end
   ρatm = atm_density
   h0 = 1e-3
 
-  ####
-  # cons2prim from Kastaun et al. arXiv:2005.01821
-  ####
+  #### cons2prim from Kastaun et al. arXiv:2005.01821
 
   γuurr = 1/γrr
   q = τ / D # (22)-(23)
@@ -186,43 +180,27 @@ end
 
     # root bracket adjustment, cf. Appendix
     if D / What_b < ρmax < D
-      # @warn "adjusting upper root bracket"
-
       # infamous closure bug ...
       # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
       fn_ρmax = let D = D, ρmax = ρmax
         mu -> D / fn_What(mu) - ρmax
       end
-
       mu_b = find_zero(fn_ρmax, (mu_a, mu_plus), Roots.A42())
-
       if isnan(mu_b)
         error("Failed to correct upper interval bound")
       end
-
-      # println("New interval: ($mu_a, $mu_b)")
     end
 
     if D / What_b < ρmin < D # yes, we need What_b here
-      # @warn "adjusting lower root bracket"
-      # @show D/What_b, ρmin, D
-
       # infamous closure bug ...
       # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
       fn_ρmin = let D = D, ρmin = ρmin
         mu -> D / fn_What(mu) - ρmin
       end
-
-      # @show(mu_a, fn_What(mu_a))
-      # @show(mu_b, fn_What(mu_b))
       mu_a = find_zero(fn_ρmin, (mu_a, mu_b), Roots.A42())
-      # @show(mu_a, fn_What(mu_a))
-
       if isnan(mu_b)
         error("Failed to correct lower interval bound")
       end
-
-      # println("New interval: ($mu_a, $mu_b)")
     end
 
     # infamous closure bug ...
@@ -259,20 +237,12 @@ end
     mf_a = master_fn(mu_a)
     mf_b = master_fn(mu_b)
     mf_ab = mf_a*mf_b
-    # if !isfinite(mf_ab) || mf_ab >= 0
-    # # if mf_ab >= 0 || !isfinite(mf_ab)
-    #   println((; D, Sr, τ, mf_a, mf_b, mu_a, mu_b, r, rcoord))
-    #   error("Master function is ill conditioned!")
-    # end
-
     mu = find_zero(master_fn, (mu_a, mu_b), Roots.A42() #= =^= algorithm 748 =#)
     if isnan(mu)
       error("Failed to find root")
     end
 
-    #####
-    # evaluate primitives
-    #####
+    #### evaluate primitives
 
     # eq (36) [Kastaun] : h*W = 1/mu
     # eq (15) [Kastaun] : ρ*W = D
@@ -297,11 +267,8 @@ end
     v  = min(v, v0)
     v2 = v^2
     vr = sign(Sr) * sqrt(v2/γuurr)
-    # vr = sign(Sr) * min(abs(mu * rr), abs(v0))
 
     # eq (68) or (40)
-    # v = min(mu * r, v0)
-    # @show v2, vr, v0
     W = 1.0 / sqrt(1.0 - v2)
     # eq (41)
     ρ = D / W
@@ -312,12 +279,6 @@ end
 
     # enforce EOS bounds
 
-    # Deppe at all also have a atmopshere treatment for this
-    # if ρatm <= ρ <= ρmin
-    #   @assert v^2 < 1e-4
-    # end
-
-    # @assert ρ > ρmin ρ
     if ρ < ρmin
       # atmosphere II
 
@@ -342,37 +303,11 @@ end
       # error()
       # adjusting τ such that energy density is within bounds
       # cf. 3rd paragraph in sec. III.A in Kastaun paper
-      # @show ϵ, ϵmin, rcoord
       ϵ = ϵmin
       # invert this for τ: ϵ = W * (q - mu * r2) + v^2 * W^2 / (1.0 + W), q = τ / D
       τ = (ϵ / W  - v2 * W / (1.0 + W) + mu * r2) * D
     end
 
-
-    # if ϵ < ϵmin
-    #   error()
-    #   # adjusting τ such that energy density is within bounds
-    #   # cf. 3rd paragraph in sec. III.A in Kasτn paper
-    #   ϵ = ϵmin
-    #   τb4 = τ
-    #   # invert this for τ: ϵ = W * (q - mu * r2) + v^2 * W^2 / (1.0 + W), q = τ / D
-    #   τ = (ϵ / W  - v2 * W / (1.0 + W) + mu * r2) * D
-    # end
-    # if ρ < ρmin || ϵ < ϵmin
-    #   @warn "atmosphering II @ r = $r : (ρ,ϵ) = ($ρ, $ϵ)"
-    #   flag_cons2prim_atm2 = 1
-    #   # atmosphering
-    #   ρ, ϵ, v, W = max(ρmin, D), ϵmin, 0.0, 1.0
-    #   D = ρ
-    #   τ =
-    #   # ρ, ϵ, v, W = ρmin, ϵmin, 0.0, 1.0
-    # elseif ρ > ρmax || ϵ > ϵmax
-    #   @warn "ρmax || ϵmax reached"
-    #   flag_cons2prim_atm2 = 1
-    #   ρ, ϵ, v = ρmax, ϵmax, Si*mu/D
-    #   D = ρ
-    # end
-
     # compute thermodynamic vars
     p = eos(Pressure, Density, ρ, InternalEnergy, ϵ)
     if p < 0