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

@@ -404,6 +404,41 @@ norm_L2(u, mesh::Mesh) = sqrt(broken_inner_product(u, u, mesh))
 norm_L1(u, mesh::Mesh) = broken_inner_product(abs.(u), ones(Float64, size(u)), mesh)
 
 
+function differentiate!(du, u, mesh::Mesh1d)
+  @unpack invdetJ = get_static_variables(mesh.cache)
+  for (duk,uk,invdetJk) in zip(eachcell(mesh,du),eachcell(mesh,u),eachcell(mesh,invdetJ))
+    differentiate!(duk, uk, mesh.element)
+    duk .*= invdetJk
+  end
+end
+function differentiate!(du, u, mesh::Mesh1d{FVElement})
+  @unpack x, invdetJ = get_static_variables(mesh.cache)
+  K, = mesh.tree.dims
+  for k in 2:K-1
+    du[k] = (u[k+1]-u[k-1])/(x[k+1]-x[k-1])
+  end
+  x1, x2, x3 = x[1], x[2], x[3]
+  M = [ 1 x1 x1^2;
+        1 x2 x2^2;
+        1 x3 x3^2 ]
+  uu = u[1:3]
+  b = M \ uu
+  du[1] = b[2] + 2*b[3]*x1
+  x1, x2, x3 = x[end], x[end-1], x[end-2]
+  M = [ 1 x1 x1^2;
+        1 x2 x2^2;
+        1 x3 x3^2 ]
+  uu = u[end:-1:end-2]
+  b = M \ uu
+  du[end] = b[2] + 2*b[3]*x1
+end
+function differentiate(u, mesh)
+  du = similar(u)
+  differentiate!(du, u, mesh)
+  return du
+end
+
+
 """
     find_nearest(xi, mesh::Mesh)
 
@@ -571,7 +606,7 @@ end
 
 struct CellDataIterator{M<:AbstractMesh,A<:AbstractArray}
   mesh::M
-  data::A
+  data::A # TODO Replace with NTuple{N,A} to allow iterating multiple arrays at once
   function CellDataIterator(m::AbstractMesh, a::AbstractArray)
     @toggled_assert length(m) == length(a)
     return new{typeof(m),typeof(a)}(m,a)