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

+#######################################################################
+#                         Conservative fixing                         #
+#######################################################################
+
+
+@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
+  ρmin = atm_density * atm_threshold
+
+  # B = 1.0 # in special relativity
+
+  #####
+  # conservative variable fixing following Deppe et al. arXiv:2109.12033
+  ######
+
+  fixed = false
+
+  if D < ρmin || τ < 0.0
+    D = ρmin
+    τ = 1e-12
+    fixed = true
+  end
+
+  # if D < 0.001 * ρmin
+  #   # @warn "fixing D $D @ r = $r"
+  #   D = ρmin
+  #   fixed = true
+  # end
+  #
+  # if τ < 0.0
+  #   # @warn "fixing τ $τ @ r = $r"
+  #   τ = 1e-12
+  #   fixed = true
+  # end
+  # fix τ such that v < 0.999
+  # if τ < 0.0
+  #   @warn "fixing τ $τ @ r = $r"
+  #   # What = 1 + τ / D = 1 / sqrt(1 - v^2)
+  #   v2 = 0.999
+  #   τ = max(D * (1.0 / sqrt(1.0 - v2) - 1), 1e-12)
+  #   display(τ)
+  #   fixed = true
+  # end
+
+  # (48)-(50) in arXiv:2109.12033 we put B=0, ̂μ = 0
+  # solution to (52) is then
+  What = 1.0 + τ / D
+
+  S2 = S^2
+  D2 = D^2
+  ϵs = 1e-12
+  factor = sqrt( (1.0-ϵs) * (What^2-1.0) * D2 / (S2 + D2 * 1e-16) )
+  # @show factor, What
+  if factor < 1.0 && What > 1.0
+    # @warn "fixing ..."
+    Sold = S
+    S = 0.99 * sign(S) * sqrt( (1.0-ϵs) * (What^2-1.0) * D2 - D2 * 1e-16 )
+    # @warn "S  = $Sold => $S @ r = $r"
+    S2 = S^2
+    factor = max(1.0, sqrt( (1.0-ϵs) * (What^2-1.0) * D2 / (S2 + D2 * 1e-16) ))
+    # @show factor
+    @assert factor >= 1.0
+    fixed = true
+  end
+
+  # h0 = 1e-3
+  # S2 = S^2 / B^2
+  # r2 = S2 / D^2
+  # z02 = r2 / h0^2
+  # v02 = z02 / (1.0 + z02)
+  # max_v2 = 0.999
+  # if v02 > max_v2
+  #   S = sign(S) * sqrt(max_v2/(1-max_v2) * h0^2 * D^2)
+  # end
+
+  flag_consfixed = Float64(fixed)
+
+  @returns D, S, τ, flag_consfixed
+end
+
+
+#######################################################################
+#                  cons2prim a la Kastaun et al 2021                  #
+#######################################################################
+
+
+@with_signature [simd=false] function cons2prim(equation::AbstractEquation)
+  @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_rescaled_spherical1d(equation::AbstractEquation)
+  @accepts D, Sr, τ, r
+  @accepts B
+  γrr = 1/B^2
+  D, Sr, τ, ρ, vr, ϵ, p = cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr)
+  # in case something was adjusted 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)
+
+  rcoord = r
+
+  @unpack atm_density, atm_threshold, eos = equation
+  ρmin = atm_density * atm_threshold
+  @unpack K, Gamma = eos
+  pmin = K * ρmin^Gamma
+  ϵmin = pmin / ((Gamma-1)*ρmin)
+  ρmax = Inf
+  ρatm = atm_density
+  h0 = 1e-3
+
+  #### cons2prim from Kastaun et al. arXiv:2005.01821
+
+  γuurr = 1/γrr
+  q = τ / D # (22)-(23)
+  rr = Sr / D # (22)-(23), Si := S_i is 3-vector
+  r2 = γuurr*rr*rr
+  r  = sqrt(r2)
+  z0 = r/h0
+  v0 = z0 / sqrt(1.0 + z0^2)
+  h02 = h0^2
+
+  if v0 > 1
+    error((v0,z0,Sr,D))
+  end
+
+  # What is not a question but rather W^hat, cf. (40)
+  fn_What = let r = r, v0 = v0
+    mu -> begin
+      vhat = min(mu * r, v0)
+      return 1.0 / sqrt(1.0 - vhat^2)
+    end
+  end
+
+  # setup bracket
+  mu_a = 0.0
+  # with no EM fields we can find the root of (49) analytically, because r = const.
+  mu_plus = 1.0 / sqrt(h02 + r2)
+  mu_b = r < h0 ? 1.0 / h0 : mu_plus
+
+  What_a = fn_What(mu_a)
+  What_b = fn_What(mu_b)
+  if D / What_b > ρmax
+    error("Invalid state encountered!")
+  end
+
+  if D < ρmin
+
+    @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
+
+  else
+
+    # root bracket adjustment, cf. Appendix
+    if D / What_b < ρmax < D
+      # infamous closure bug ...
+      # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
+      fn_ρmax = let D = D, ρmax = ρmax
+        mu -> D / fn_What(mu) - ρmax
+      end
+      mu_b = find_zero(fn_ρmax, (mu_a, mu_plus), Roots.A42())
+      if isnan(mu_b)
+        error("Failed to correct upper interval bound")
+      end
+    end
+
+    if D / What_b < ρmin < D # yes, we need What_b here
+      # infamous closure bug ...
+      # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
+      fn_ρmin = let D = D, ρmin = ρmin
+        mu -> D / fn_What(mu) - ρmin
+      end
+      mu_a = find_zero(fn_ρmin, (mu_a, mu_b), Roots.A42())
+      if isnan(mu_b)
+        error("Failed to correct lower interval bound")
+      end
+    end
+
+    # infamous closure bug ...
+    # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
+    master_fn = let D = D, v0 = v0, r = r, q = q, r2 = r2
+      mu -> begin # eq (44)
+
+        _vhat = min(mu * r, v0)
+        _What = 1.0 / sqrt(1.0 - _vhat^2)
+        ρhat = D / _What
+
+        # eq (42)
+        ϵhat = _What * (q - mu * r2) + _vhat^2 * _What^2 / (1.0 + _What)
+
+        # TODO enforce EOS bounds - needed?
+
+        # compute thermodynamic vars
+        phat = eos(Pressure, Density, ρhat, InternalEnergy, ϵhat)
+        # hhat = eos(Enthalpy, Density, ρhat, InternalEnergy, ϵhat)
+        ahat = phat / (ρhat * (1.0 + ϵhat))
+
+        # eq (46)-(48)
+        nu_A = (1.0 + ahat) * (1.0 + ϵhat) / _What
+        nu_B = (1.0 + ahat) * (1.0 + q - mu * r2)
+        nu = max(nu_A, nu_B)
+
+        # eq (44)
+        f = mu - 1.0 / (nu + r2 * mu)
+        return f
+      end
+    end
+
+    # verify if there is a solution inside bracket
+    mf_a = master_fn(mu_a)
+    mf_b = master_fn(mu_b)
+    mf_ab = mf_a*mf_b
+    mu = find_zero(master_fn, (mu_a, mu_b), Roots.A42() #= =^= algorithm 748 =#)
+    if isnan(mu)
+      error("Failed to find root")
+    end
+
+    #### evaluate primitives
+
+    # eq (36) [Kastaun] : h*W = 1/mu
+    # eq (15) [Kastaun] : ρ*W = D
+    # ρhW2 = D / mu
+
+    # regularize velocity according to Dumbser et al arxiv:2307.06629
+    # eqs (82)-(84)
+    # y, y0 = ρhW2, 1e-4
+    # if y > y0
+      # v = S / mu
+      # compute v accoriding to eq (68) or (40) [Kasτn]
+      # vr = sign(Sr) * min(mu * r, v0)
+      # vr = sign(Sr) * min(abs(mu * rr), abs(v0))
+    # else
+    #   # @warn "limiting v"
+    #   # limit velocity ratio v/S
+    #   fy = 2.0 * (y/y0)^3 - 3.0 * (y/y0)^2 + 1.0
+    #   v = Si * y / (y^2 + fy * 5e-9)
+    # end
+
+    v = sqrt(mu^2 * γuurr*rr*rr)
+    v  = min(v, v0)
+    v2 = v^2
+    vr = sign(Sr) * sqrt(v2/γuurr)
+
+    # eq (68) or (40)
+    W = 1.0 / sqrt(1.0 - v2)
+    # eq (41)
+    ρ = D / W
+    # eq (42)
+    ϵ = W * (q - mu * r2) + v2 * W^2 / (1.0 + W)
+    # eq (38)
+    # ϵ2 = W * (q - mu * r2) + W - 1
+
+    # enforce EOS bounds
+
+    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 Kastaun paper
+      ϵ = ϵmin
+      # 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
+
+    # compute thermodynamic vars
+    p = eos(Pressure, Density, ρ, InternalEnergy, ϵ)
+    if p < 0
+      error((ρ,ϵ))
+    end
+
+  end
+
+  return D, Sr, τ, ρ, vr, ϵ, p
+end