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

@@ -154,6 +154,145 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spheric
 end
 
 
+function initialdata_equation!(::Val{:tov}, env, P::Project{:spherical1d}, prms)
+
+  filename = prms["id_filename"]
+  if !isfile(filename)
+    error("GRHD: Failed to locate TOV initial data at id_filename = $filename")
+  end
+
+  data = h5open(filename, "r") do file
+    g = file["TOV"]
+    Γ  = read_attribute(g, "Γ")
+    K  = read_attribute(g, "K")
+    ρc = read_attribute(g, "ρc")
+    M  = read_attribute(g, "M")
+    rsch = read(g, "r") # schwarzschild radius coordinate
+    riso = read(g, "R") # isotropic radius coordinate
+    ρ   = read(g, "ρ")
+    p   = read(g, "p")
+    ϵ   = read(g, "ϵ")
+    ϕ   = read(g, "ϕ")
+    m   = read(g, "m")
+
+    # metric potentials: ds^2 = -A^2 dt^2 + B^2 γij dx^i dx^j
+    A = @. exp(ϕ)
+    B = @. rsch / riso
+    # B[1] is nan, because rsch = riso = 0 there
+    # use a parabola to interpolate data point at riso = 0
+    # this also ensures that B is symmetric arcoss riso = 0
+    B2, B3, r2, r3 = B[2], B[3], riso[2], riso[3]
+    B[1] = (r3^2*B2-r2^2*B3)/(r3^2-r2^2)
+
+    (; Γ, K, ρc, M, riso, ρ, ϵ, p, A, B)
+  end
+
+  @unpack equation = P
+  @unpack mesh = env
+  @unpack eos = P.equation
+  @assert data.K ≈ eos.K
+  @assert data.Γ ≈ eos.Gamma
+
+  equation.atm_density = data.ρc * equation.atm_factor
+  equation.atm_factor = prms["atm_factor"]
+  equation.atm_threshold = prms["atm_threshold_factor"]
+
+  @unpack D, Sr, τ = get_dynamic_variables(mesh.cache)
+  @unpack ρ, p, ϵ, vr, r, max_v,
+          init_D, init_Sr, init_τ, init_vr, init_p,
+          init_max_v = get_static_variables(mesh.cache)
+  @unpack α, βru, ∂αr, ∂βrrdu,
+          γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
+          Γrrr, Γrθθ, Γrϕϕ, Γθrθ, Γθϕϕ, Γϕrϕ, Γϕθϕ,
+          Krr, Kθθ, Kϕϕ = get_static_variables(mesh.cache)
+
+  if minimum(r) ≈ 0.0
+    error("gird includes coordinate origin, but spherical1d formulation is singular at r=0")
+  end
+
+  # need to split initial data into portion that lies inside and outside of star,
+  # because derivatives of the hydro variables are discontinuous at the surface
+  # although metric functions are supposed to be better conditioned there,
+  # we also split them for the interpolation, because of potential contamination
+  # with numerical error from patching
+  rmax = data.riso[end]
+  idx_exterior = findfirst(ri -> ri > rmax, r) |> something
+  @assert idx_exterior > 0
+  interior = 1:idx_exterior-1
+  exterior = idx_exterior:length(mesh)
+
+  rint = view(r, interior)
+  rext = view(r, exterior)
+  A, B, ∂Ar, ∂Br = [ similar(r) for _ in 1:4 ]
+
+  ρ[interior] .= interpolate(rint, data.riso, data.ρ)
+  ϵ[interior] .= interpolate(rint, data.riso, data.ϵ)
+  p[interior] .= interpolate(rint, data.riso, data.p)
+  A[interior] .= interpolate(rint, data.riso, data.A)
+  B[interior] .= interpolate(rint, data.riso, data.B)
+
+  ρ[exterior] .= equation.atm_density
+  p[exterior] .= @. eos(Pressure, Density, ρ[exterior])
+  ϵ[exterior] .= @. eos(InternalEnergy, Density, ρ[exterior])
+  # metric potentials of Schwarzschild solution in isotropic coordinates:
+  #   ds^2 = -A^2 dt^2 + B^2 γij dx^i dx^j
+  @unpack M = data
+  A[exterior] .= @. (1-M/(2*rext))/(1+M/(2*rext))
+  B[exterior] .= @. (1+M/(2*rext))^2
+
+  # BasicInterpolators.jl cannot differentiate ...  so we diff numerically
+  ∂Ar .= dg1d.differentiate(A, mesh)
+  ∂Br .= dg1d.differentiate(B, mesh)
+
+  θ = π/2
+
+  @. α = A
+  @. βru = 0.0
+  @. ∂αr = ∂Ar
+  @. ∂βrrdu = 0
+
+  @. γrr = B^2
+  @. γθθ = B^2*r^2
+  @. γϕϕ = B^2*r^2*sin(θ)^2
+
+  @. ∂γrrr = 2*B*∂Br
+  @. ∂γrθθ = 2*(r*B^2+B*∂Br*r^2)
+  @. ∂γrϕϕ = ∂γrθθ * sin(θ)^2
+  @. ∂γθϕϕ = 0.0
+
+  @. Γrrr = ∂Br/B
+  @. Γrθθ = (-r^2*B*∂Br-r*B^2)/B^2
+  @. Γrϕϕ = Γrθθ*sin(θ)^2
+  @. Γθrθ = (r^2*B*∂Br+r*B^2)/(r^2*B^2)
+  @. Γθϕϕ = -sin(θ)*cos(θ)
+  @. Γϕrϕ = (r^2*B*∂Br+r*B^2)/(r^2*B^2)
+  @. Γϕθϕ = cos(θ)/sin(θ)
+
+  @. Krr = 0.0
+  @. Kθθ = 0.0
+  @. Kϕϕ = 0.0
+
+  @. vr = 0
+  W = 1.0
+  h = @. 1 + ϵ + p/ρ
+  @. D  = W * ρ
+  @. Sr = W^2 * ρ * h * vr
+  @. τ  = W^2 * ρ * h - p - D
+
+  broadcast_volume!(cons2prim, equation, mesh)
+  broadcast_volume!(maxspeed, equation, mesh)
+  broadcast_volume!(flux_source_spherical1d, equation, mesh)
+
+  @. init_D     = D
+  @. init_Sr    = Sr
+  @. init_τ     = τ
+  @. init_vr    = vr
+  @. init_p     = p
+  @. init_max_v = max_v
+
+end
+
+
 function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1d}, prms)
 
   filename = prms["id_filename"]