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

@@ -421,6 +421,68 @@ function impose_symmetry!(equations::Equation, mesh::Mesh)
   return
 end
 
+function impose_symmetry!(equations::Equation, mesh::Mesh{FVElement})
+
+  @unpack cache = mesh
+  @unpack rD = get_dynamic_variables(cache)
+  @unpack D, Sr, τ, r = get_static_variables(cache)
+  L = layout(mesh)
+  Npts, _ = L
+
+  @show rD[1], rD[2]
+
+  # symmetric variables are interpolated such that derivative vanishes at r=0
+  for var in (D, τ)
+    tmp = 0
+    # use quadratic interpolation around origin
+    var2, var3, r2, r3 = var[2], var[3], r[2], r[3]
+    var[1] = (r3^2*var2-r2^2*var3)/(r3^2-r2^2)
+  end
+  # asymmetric variables are set to zero at r=0
+  Sr[1] = 0
+
+  @show D[1], D[2]
+
+  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
+  @accepts A, B, ∂Br, r
+  @accepts [bdry] nx
+
+  γuurr    = 1/B^2
+  init_vru = γuurr * init_vr
+  vru      = γuurr * vr
+  B3       = B^3
+
+  # # Dirichlet conditions
+  bdry_D   = -D   + 2*init_D
+  bdry_Sr  = -Sr  + 2*init_Sr
+  bdry_τ   = -τ   + 2*init_τ
+  bdry_vru = -vru + 2*init_vru
+  bdry_p   = -p   + 2*init_p
+  # bdry_D   = init_D
+  # bdry_Sr  = init_Sr
+  # bdry_τ   = init_τ
+  # bdry_vru = init_vru
+  # bdry_p   = init_p
+
+  @show init_D, D
+
+  nflr_rD  = r*A*B3* bdry_D * bdry_vru
+  nflr_rSr = r*A*B3* (bdry_Sr * bdry_vru + bdry_p)
+  nflr_rτ  = r*A*B3* (bdry_τ + bdry_p) * bdry_vru
+
+  # FIXME We are abusing the bdry_D,Sr,τ buffers to store the rescaled bdry values
+  bdry_D  *= r*B3
+  bdry_Sr *= r*B3
+  bdry_τ  *= r*B3
+
+  @returns [bdry] bdry_D, bdry_Sr, bdry_τ, nflr_rD, nflr_rSr, nflr_rτ
+end
+
 # @with_signature function flux(equation::Equation)
 #   @accepts sD, sSr, stau, p, rho, vr, eps, A, B, r
 #   rAB3 = r * A * B^3