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

@@ -464,6 +464,28 @@ function impose_symmetry!(P::Project{:rescaled_spherical1d}, mesh::Mesh{FVElemen
   return
 end
 
+
+function impose_symmetry_sources!(P::Project{:spherical1d}, mesh::Mesh)
+  @unpack r = get_static_variables(mesh.cache)
+  r[1] ≈ 0.0 || return
+  @unpack src_D, src_Sr, src_τ = get_static_variables(mesh.cache)
+  L = layout(mesh)
+  Npts, _ = L
+  v_src_D  = dg1d.vreshape(src_D, L)
+  v_src_τ  = dg1d.vreshape(src_τ, L)
+  # even variables are interpolated such that derivative vanishes at r=0
+  for var in (v_src_D, v_src_τ)
+    tmp = 0
+    for j = 2:Npts
+      tmp += mesh.element.D[1,j] * var[j,1]
+    end
+    var[1,1] = -tmp/mesh.element.D[1,1]
+  end
+  # odd variables are set to zero at r=0
+  src_Sr[1] = 0
+  return
+end
+
 @with_signature function fv_bdry_flux_rescaled_spherical1d(equation::Equation)
   @accepts D, Sr, τ, vr, p
   @accepts init_D, init_Sr, init_τ, init_vr, init_p