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

@@ -6,6 +6,8 @@
 @with_signature function conservative_fixing(equation::AbstractEquation)
   @accepts D, S, τ, r
 
+  TODO()
+
   @unpack atm_density, atm_threshold, eos = equation
   # atm_density = equation.atmosphere.ρ
   # atm_threshold = equation.atmosphere.threshold
@@ -89,8 +91,29 @@ end
 
 
 @with_signature [simd=false] function cons2prim_kastaun(equation::AbstractEquation)
-  @accepts D, Sr, τ
+  @accepts D, Sr, τ, r
   @accepts γrr
+  D, Sr, τ, ρ, vr, ϵ, p = cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr)
+  @returns D, Sr, τ, ρ, vr, ϵ, p
+end
+@with_signature [simd=false] function cons2prim_kastaun_rescaled(equation::AbstractEquation)
+  @accepts D, Sr, τ, r
+  @accepts A, B
+  γrr = 1/B^2
+  D, Sr, τ, ρ, vr, ϵ, p = cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr)
+  # in case something was adjustbed by cons2prim
+  rD  = r*B^3*D
+  rSr = r*B^3*Sr
+  rτ  = r*B^3*τ
+  @returns D, Sr, τ, rD, rSr, rτ, ρ, vr, ϵ, p
+end
+
+
+@inline function cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr)
+  # @accepts D, Sr, τ, r
+  # @accepts γrr
+
+  rcoord = r
 
   @unpack atm_density, atm_threshold, eos = equation
   # atm_density = equation.atmosphere.ρ
@@ -140,12 +163,9 @@ end
   end
 
   if D < ρmin
-    @warn "atmosphering I"
 
-    error((D,Sr,τ))
-
-    @assert eos isa EquationOfState.IdealGas
-    @unpack Gamma = eos
+    @assert eos isa EquationOfState.Polytrope
+    @unpack K, Gamma = eos
 
     # employ atmosphere values by taking the zero temperature limit
     # for an ideal gas eos, we have that eos.Gamma becomes the Gamma for
@@ -154,25 +174,18 @@ end
     ρ = ρatm
     vr = 0.0
     W = 1.0
-    # p = K * ρ^Gamma
-    # ϵ = Gamma * ρ^
-    # D = W * ρ
-    # S = ρ * (1+ϵ+p/ρ) * W^2 * v
-    # τ = ρ * (1+ϵ+p/ρ) * W^2 - p - D
-
-    # if eos isa ColdEquationOfState
-    #   ϵ = eos.cold_eos(InternalEnergy, Density, ρ)
-    #   p = eos.cold_eos(Pressure, Density, ρ)
-    # else
-    #   ϵ = eos(InternalEnergy, Density, ρ)
-    #   p = eos(Pressure, Density, ρ)
-    # end
+    p = K * ρ^Gamma
+    ϵ = p / ((Gamma-1)*ρ)
+    h = 1 + ϵ + p/ρ
+    D  = W*ρ
+    Sr = W^2*ρ*h*vr
+    τ  = W^2*ρ*h - p - D
 
   else
 
     # root bracket adjustment, cf. Appendix
     if D / What_b < ρmax < D
-      @warn "adjusting upper root bracket"
+      # @warn "adjusting upper root bracket"
 
       # infamous closure bug ...
       # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
@@ -186,11 +199,12 @@ end
         error("Failed to correct upper interval bound")
       end
 
-      println("New interval: ($mu_a, $mu_b)")
+      # println("New interval: ($mu_a, $mu_b)")
     end
 
     if D / What_b < ρmin < D # yes, we need What_b here
-      @warn "adjusting lower root bracket"
+      # @warn "adjusting lower root bracket"
+      # @show D/What_b, ρmin, D
 
       # infamous closure bug ...
       # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
@@ -198,13 +212,16 @@ end
         mu -> D / fn_What(mu) - ρmin
       end
 
+      # @show(mu_a, fn_What(mu_a))
+      # @show(mu_b, fn_What(mu_b))
       mu_a = find_zero(fn_ρmin, (mu_a, mu_b), Roots.A42())
+      # @show(mu_a, fn_What(mu_a))
 
       if isnan(mu_b)
         error("Failed to correct lower interval bound")
       end
 
-      println("New interval: ($mu_a, $mu_b)")
+      # println("New interval: ($mu_a, $mu_b)")
     end
 
     # infamous closure bug ...
@@ -241,10 +258,11 @@ end
     mf_a = master_fn(mu_a)
     mf_b = master_fn(mu_b)
     mf_ab = mf_a*mf_b
-    if mf_a * mf_b >= 0 || !isfinite(mf_ab)
-      println((; D, Sr, τ, mf_a, mf_b, mu_a, mu_b, r))
-      error("Master function is ill conditioned!")
-    end
+    # if !isfinite(mf_ab) || mf_ab >= 0
+    # # if mf_ab >= 0 || !isfinite(mf_ab)
+    #   println((; D, Sr, τ, mf_a, mf_b, mu_a, mu_b, r, rcoord))
+    #   error("Master function is ill conditioned!")
+    # end
 
     mu = find_zero(master_fn, (mu_a, mu_b), Roots.A42() #= =^= algorithm 748 =#)
     if isnan(mu)
@@ -266,7 +284,7 @@ end
       # v = S / mu
       # compute v accoriding to eq (68) or (40) [Kasτn]
       # vr = sign(Sr) * min(mu * r, v0)
-      vr = mu * rr
+      # vr = sign(Sr) * min(abs(mu * rr), abs(v0))
     # else
     #   # @warn "limiting v"
     #   # limit velocity ratio v/S
@@ -274,9 +292,15 @@ end
     #   v = Si * y / (y^2 + fy * 5e-9)
     # end
 
+    v = sqrt(mu^2 * γuurr*rr*rr)
+    v  = min(v, v0)
+    vr = sign(Sr) * sqrt(v^2/γuurr)
+    # vr = sign(Sr) * min(abs(mu * rr), abs(v0))
+
     # eq (68) or (40)
     # v = min(mu * r, v0)
-    v2 = γuurr*vr*vr
+    v2 = v^2
+    # @show v2, vr, v0
     W = 1.0 / sqrt(1.0 - v2)
     # eq (41)
     ρ = D / W
@@ -292,9 +316,29 @@ end
     #   @assert v^2 < 1e-4
     # end
 
-    @assert ρ > ρmin
-    if ϵ < ϵmin
-      error()
+    # @assert ρ > ρmin ρ
+    if ρ < ρmin
+      # atmosphere II
+
+      @assert eos isa EquationOfState.Polytrope
+      @unpack K, Gamma = eos
+
+      # employ atmosphere values by taking the zero temperature limit
+      # for an ideal gas eos, we have that eos.Gamma becomes the Gamma for
+      # a polytrope, see grhd notes for derivation
+
+      ρ = ρatm
+      vr = 0.0
+      W = 1.0
+      p = K * ρ^Gamma
+      ϵ = p / ((Gamma-1)*ρ)
+      h = 1 + ϵ + p/ρ
+      D  = W*ρ
+      Sr = W^2*ρ*h*vr
+      τ  = W^2*ρ*h - p - D
+
+    elseif ϵ < ϵmin
+      # error()
       # adjusting τ such that energy density is within bounds
       # cf. 3rd paragraph in sec. III.A in Kasτn paper
       ϵ = ϵmin
@@ -302,6 +346,17 @@ end
       # invert this for τ: ϵ = W * (q - mu * r2) + v^2 * W^2 / (1.0 + W), q = τ / D
       τ = (ϵ / W  - v2 * W / (1.0 + W) + mu * r2) * D
     end
+
+
+    # if ϵ < ϵmin
+    #   error()
+    #   # adjusting τ such that energy density is within bounds
+    #   # cf. 3rd paragraph in sec. III.A in Kasτn paper
+    #   ϵ = ϵmin
+    #   τb4 = τ
+    #   # invert this for τ: ϵ = W * (q - mu * r2) + v^2 * W^2 / (1.0 + W), q = τ / D
+    #   τ = (ϵ / W  - v2 * W / (1.0 + W) + mu * r2) * D
+    # end
     # if ρ < ρmin || ϵ < ϵmin
     #   @warn "atmosphering II @ r = $r : (ρ,ϵ) = ($ρ, $ϵ)"
     #   flag_cons2prim_atm2 = 1
@@ -325,6 +380,5 @@ end
 
   end
 
-  # @returns D, S, τ, ρ, v, ϵ, p
-  @returns τ, ρ, vr, ϵ, p
+  return D, Sr, τ, ρ, vr, ϵ, p
 end