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/equation.jl

@@ -132,12 +132,30 @@ end
 end
 
 
+# @with_signature function maxspeed(equation::Equation)
+#   @accepts vr
+#   @accepts α, βru, γrr
+#   γuurr = 1/γrr
+#   vru = γuurr * vr
+#   max_v = abs(α*vru - βru)
+#   @returns max_v
+# end
 @with_signature function maxspeed(equation::Equation)
-  @accepts vr
-  @accepts α, βru, γrr
+  @accepts vr, ρ, ϵ
+  @accepts α, γrr
+  @unpack eos = equation
+  cs2 = eos(SquareSpeedOfSound, Density, ρ, InternalEnergy, ϵ)
+  cs = sqrt(abs(cs2))
   γuurr = 1/γrr
   vru = γuurr * vr
-  max_v = abs(α*vru - βru)
+  v2 = vru*vr
+  A = α
+  λ0 = A*vru
+  W = 1/sqrt(1-v2)
+  λp = A/(1-v2)*(vru*(1-cs2)+cs/W*sqrt((1-v2*cs2)*γuurr-vru^2*(1-cs2)))
+  λm = A/(1-v2)*(vru*(1-cs2)-cs/W*sqrt((1-v2*cs2)*γuurr-vru^2*(1-cs2)))
+  max_v = absmax(λm, λ0, λp)
+  # TODO Can't use this method with bondi accretion test, need to use old version ...
   @returns max_v
 end