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

@@ -125,6 +125,47 @@ function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
 end
 
 
+function rhs!(mesh::Mesh1d, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
+
+  @unpack cache = mesh
+  @unpack equation = P
+
+  @unpack rD, rSr, rτ                    = get_dynamic_variables(cache)
+  @unpack rhs_rD, rhs_rSr, rhs_rτ        = get_rhs_variables(cache)
+  @unpack flr_rD, flr_rSr, flr_rτ,
+          max_v, vr, p, D, Sr, τ,
+          src_rD, src_rSr, src_rτ        = get_static_variables(cache)
+  @unpack nflr_rD, nflr_rSr, nflr_rτ,
+          bdry_D, bdry_Sr, bdry_τ,
+          bdry_max_v, bdry_vr, bdry_p = get_bdry_variables(cache)
+
+  # @show(extrema(rD))
+  # @show(extrema(rSr))
+  # @show(extrema(rτ))
+  broadcast_volume!(unscale_conservatives, equation, mesh)
+  impose_symmetry!(equation, mesh)
+  broadcast_volume!(cons2prim_kastaun_rescaled, equation, mesh)
+  broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, D,     bdry_D)
+  dg1d.interpolate_face_data!(mesh, Sr,    bdry_Sr)
+  dg1d.interpolate_face_data!(mesh, τ,     bdry_τ)
+  dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
+  dg1d.interpolate_face_data!(mesh, vr,    bdry_vr)
+  dg1d.interpolate_face_data!(mesh, p,     bdry_p)
+
+  broadcast_volume!(flux_source_rescaled_spherical1d, equation, mesh)
+  broadcast_faces!(llf_rescaled_spherical1d, equation, mesh)
+  broadcast_bdry!(bdryllf_rescaled_spherical1d, equation, mesh)
+
+  compute_rhs_weak_form!(rhs_rD,  flr_rD,  src_rD,  nflr_rD,  mesh)
+  compute_rhs_weak_form!(rhs_rSr, flr_rSr, src_rSr, nflr_rSr, mesh)
+  compute_rhs_weak_form!(rhs_rτ,  flr_rτ,  src_rτ,  nflr_rτ,  mesh)
+
+  return
+end
+
+
 function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction)
 
   TODO()