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

@@ -98,10 +98,10 @@ end
 end
 @with_signature [simd=false] function cons2prim_rescaled_spherical1d(equation::AbstractEquation)
   @accepts D, Sr, τ, r
-  @accepts A, B
+  @accepts B
   γrr = 1/B^2
   D, Sr, τ, ρ, vr, ϵ, p = cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr)
-  # in case something was adjustbed by cons2prim
+  # in case something was adjusted by cons2prim
   rD  = r*B^3*D
   rSr = r*B^3*Sr
   rτ  = r*B^3*τ
@@ -294,12 +294,12 @@ end
 
     v = sqrt(mu^2 * γuurr*rr*rr)
     v  = min(v, v0)
-    vr = sign(Sr) * sqrt(v^2/γuurr)
+    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)
-    v2 = v^2
     # @show v2, vr, v0
     W = 1.0 / sqrt(1.0 - v2)
     # eq (41)
@@ -340,9 +340,8 @@ end
     elseif ϵ < ϵmin
       # error()
       # adjusting τ such that energy density is within bounds
-      # cf. 3rd paragraph in sec. III.A in Kasτn paper
+      # cf. 3rd paragraph in sec. III.A in Kastaun 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