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

math/3plus1.py

@@ -7,7 +7,7 @@
 #       extension: .py
 #       format_name: light
 #       format_version: '1.5'
-#       jupytext_version: 1.12.0
+#       jupytext_version: 1.16.0
 #   kernelspec:
 #     display_name: Python 3 (ipykernel)
 #     language: python
@@ -95,6 +95,10 @@ display(γ4ud)
 Γs = christoffels(γ_coords, γ, γinv)
 display(Γs)
 
+# $$
+#     \gamma_{ij} = \eta_{ij}
+# $$
+
 η_coords = coords[1:]
 display(η_coords)
 η = sp.diag(1, r**2, r**2 * sp.sin(θ)**2)
@@ -147,6 +151,20 @@ display(γhat.det())
 γhat_Γs = christoffels(γhat_coords, γhat, γhatinv)
 display(γhat_Γs)
 
+# $$
+#     \gamma_{ij} = r^{-8/3} \eta_{ij}
+# $$
+
+γhat_coords = coords[1:]
+display(γhat_coords)
+γhat = r**(sp.Rational(-8,3)) * η
+γhatinv = γhat.inv()
+display(γhat)
+display(γhatinv)
+display(γhat.det())
+γhat_Γs = christoffels(γhat_coords, γhat, γhatinv)
+display(γhat_Γs)
+
 
 # ## Extrinsic curvature
 #
@@ -195,7 +213,7 @@ for (c, Γ) in zip(sphere_coords, sphere_Γs):
             if Γij:
                 ci = sphere_coords[i];
                 cj = sphere_coords[j];
-                display(Math(f"\Gamma^{{ {c} }}_{{ {ci} {cj} }} = " + sp.latex(Γij)))
+                display(Math(fr"\Gamma^{{ {c} }}_{{ {ci} {cj} }} = " + sp.latex(Γij)))
                 # print(sp.latex(Γij))