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

@@ -182,9 +182,9 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
     # use a parabola to interpolate data point at riso = 0
     # this also ensures that B is symmetric arcoss riso = 0
     B2, B3, r2, r3 = B[2], B[3], riso[2], riso[3]
-    B[1] = B2 - r2^2*(B3-B2)/(B3^2-B2^2)
+    B[1] = (r3^2*B2-r2^2*B3)/(r3^2-r2^2)
 
-    (; Γ, K, ρc, M, rsch, riso, ρ, ϵ, p, A, B)
+    (; Γ, K, ρc, M, riso, ρ, ϵ, p, A, B)
   end
 
   @unpack equation = P
@@ -207,22 +207,16 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
   # need to split initial data into portion that lies inside and outside of star,
   # because derivatives of the hydro variables are discontinuous at the surface
   # although metric functions are supposed to be better conditioned there,
-  # we also split them for the interpolation, because of potential contamination of numerical error
+  # we also split them for the interpolation, because of potential contamination
+  # with numerical error from patching
   rmax = data.riso[end]
   idx_exterior = findfirst(ri -> ri > rmax, r) |> something
   @assert idx_exterior > 0
   interior = 1:idx_exterior-1
   exterior = idx_exterior:length(mesh)
 
-  @show r[idx_exterior-1], r[idx_exterior]
-
   rint = view(r, interior)
   rext = view(r, exterior)
-  rsch = similar(r)
-  # rsch, A, ∂Ar, B, ∂Br = [ similar(r) for _ in 1:5 ]
-
-  rsch[interior] .= interpolate(rint, data.riso, data.rsch)
-  @. rsch[exterior] = 1/2*(rext - M + sqrt(rext^2 - 2*M*rext))
 
   ρ[interior] .= interpolate(rint, data.riso, data.ρ)
   ϵ[interior] .= interpolate(rint, data.riso, data.ϵ)
@@ -230,13 +224,9 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
   A[interior] .= interpolate(rint, data.riso, data.A)
   B[interior] .= interpolate(rint, data.riso, data.B)
 
-  @show equation.atm_density
   ρ[exterior] .= equation.atm_density
   p[exterior] .= @. eos(Pressure, Density, ρ[exterior])
   ϵ[exterior] .= @. eos(InternalEnergy, Density, ρ[exterior])
-  @show extrema(ρ[exterior])
-  @show extrema(p[exterior])
-  @show extrema(ϵ[exterior])
   # metric potentials of Schwarzschild solution in isotropic coordinates:
   #   ds^2 = -A^2 dt^2 + B^2 γij dx^i dx^j
   @unpack M = data
@@ -253,11 +243,9 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
   @. D  = W * ρ
   @. Sr = W^2 * ρ * h * vr
   @. τ  = W^2 * ρ * h - p - D
-  @. rD  = r*A*B^3*D
-  @. rSr = r*A*B^3*Sr
-  @. rτ  = r*A*B^3*τ
-
-  @show(extrema(D))
+  @. rD  = r*B^3*D
+  @. rSr = r*B^3*Sr
+  @. rτ  = r*B^3*τ
 
   broadcast_volume!(unscale_conservatives, equation, mesh)
   impose_symmetry!(equation, mesh)