Skip to content
Snippets Groups Projects

Refactor GRHD project

Merged Refactor GRHD project
fa/grhd into main
Merged Florian Atteneder requested to merge fa/grhd into main
Compare
and
Show latest version
6 files
+ 90
− 891
Compare changes
  • Side-by-side
  • Inline
Files
6
src/GRHD/cons2prim.old.jl deleted 100644 → 0
+ 0
− 320
@with_signature function cons2prim(equation::Equation{<:HotEquationOfState})
@accepts D, Sr, τ, p, r, B
error()
# S2 = (Sr/B)^2
# @unpack eos, cons2prim_newtonsolver, cons2prim_bisection = equation
# p, converged = find_pressure_insane(D, S2, τ, p, r, eos,
# cons2prim_newtonsolver, cons2prim_bisection)
# if !converged && eltype(equation) == HybridEoS
# rho, success = find_density_insane(D, S2, p, r, equation, eos,
# cons2prim_newtonsolver, cons2prim_bisection)
# if !success
# error()
# end
# p = eos(Pressure, Density, rho)
# vr = Formulae.vr(B, D, Sr, τ, p)
# ϵ = Formulae.ϵ(B, D, Sr, τ, p)
# else
# rho = Formulae.ρ(B, D, Sr, τ, p)
# vr = Formulae.vr(B, D, Sr, τ, p)
# ϵ = Formulae.ϵ(B, D, Sr, τ, p)
# end
D, Sr, τ = conservative_fixing(D, Sr, τ, B, equation)
rho, vr, ϵ, p = cons2prim_kasτn((D, Sr, τ, p, r, B), equation)
# rho, vr, ϵ, p, D, Sr, τ = atmosphering(rho, vr, ϵ, p, D, Sr, τ)
@returns rho, vr, ϵ, p, D, Sr, τ
end
@with_signature function cons2prim_kastaun(equation::Equation{<:ColdEquationOfState})
@accepts D, Sr, τ, p, r
S2 = (Sr/B)^2
rho, success = find_density_insane(D, S2, p, equation)
@unpack eos = equation
p = eos(Pressure, Density, rho)
vr = Formulae.vr(B, D, Sr, τ, p)
ϵ = Formulae.ϵ(B, D, Sr, τ, p)
@returns rho, vr, ϵ, p
end
#######################################################################
# Insane version - bail out if anything is wrong #
#######################################################################
function find_pressure_insane(D, S2, τ, p, r, eos::HotEquationOfState,
cons2prim_newtonsolver, cons2prim_bisection)
# upper limit on particle velocity
# omv2 = 1e-3 # = 1 - v^2
# pmin1 = - (τ + D) + sqrt(S2 + omv2) # √((τ + D + p)^2 - S^2) >= √(1 - vmax^2)
pmin1 = - (τ + D) + sqrt(S2) + 1e-10 # √((τ + D + p)^2 - S^2) >= √(1 - vmax^2)
pmin2 = 1e-10 # physical constraint
pmin = max(pmin1, pmin2)
floor_pressure(pnew) = max(pnew, pmin)
# p0 = 0.5 * (eos.gamma - 1) * (τ + D)
p0 = p
f = pp -> pressure_equation_insane(D, S2, τ, pp, r, eos)
p, err, converged, tol = find_root(f, p0, cons2prim_newtonsolver, sanitize=floor_pressure)
if !converged
# reset to guess
p = p0
end
p = max(pmin, p0)
return p, converged
end
function pressure_equation_insane(D, S2, τ, p, r, eos)
a = τ + p + D
a2mS2 = a^2 - S2
if a2mS2 < 0
return NaN
end
W = a / sqrt(a2mS2)
rho = D / W
# if rho < 0 #1e-15
# # return NaN
# # rho = 1e-10
# end
ϵ_rho = (a - W^2 * p) / W^2 - rho
if ϵ_rho < 0
# in low temperature limit we should switch to one parameter EoS
return NaN
end
ϵ = ϵ_rho / rho
eos_p = eos(Pressure, Density, rho, InternalEnergy, ϵ) + 1e-10
f = eos_p - p
return f
end
function find_density_insane(D, S2, rho, eos,
cons2prim_newtonsolver, cons2prim_bisection)
rhomin = 1e-10
rho_init = rho
floor_density(rhonew) = max(rhonew, rhomin)
g = rho -> density_equation(D, S2, rho, eos)
rho, err, converged, tol = find_root(g, rho_init, cons2prim_newtonsolver,
sanitize=floor_density)
if !converged
# try bisection
# Figure out brackets for a bisection method.
# if bisection failed reset
rho = rho_init
end
success = converged
return rho, success
end
function density_equation(D, S2, rho, eos)
D2 = D^2
h = eos(Enthalpy, Density, rho)
W = sqrt(1 + S2 / (D2 * h^2))
chi = eos(Derivative{Pressure,Density}, Density, rho)
g = W * rho - D
# assuming zero temperature limit: dϵ/drho = P / rho^2
dg = W - S2 * chi / (W * D2 * h^3)
return g, dg
end
#######################################################################
# cons2prim a la Kasτn et al 2021 #
#######################################################################
@with_signature function cons2prim_kasτn(equation::AbstractEquation)
@accepts D, Sr, τ, p, r, B
@unpack eos, atm_density, atm_threshold = equation
rhomin, rhomax = atm_density * atm_threshold, Inf
fn_ϵmin, fn_ϵmax = rho -> -(1-1e-10), rho -> Inf
h0 = 1e-3
rhoatm = atm_density
q = τ / D
S2 = (Sr/B)^2
s = sqrt(S2) / D # =^= r in paper, r already denotes radius here
z0 = s / h0
v0 = z0 / sqrt(1 + z0^2)
h02 = h0^2
s2 = s^2
# What is not a question but rather W^hat, cf. (40)
function fn_What(mu)
_vhat = min(mu * s, v0)
if _vhat < 0
error()
println((; q, _vhat, s, mu, v0))
end
return 1 / sqrt(1 - _vhat^2)
end
# setup bracket
mu_a = 0.0
# with no EM fields we can find the root of (49) analytically, because s = const.
mu_plus = 1 / sqrt(h02 + s2)
mu_b = r < h0 ? 1 / h0 : mu_plus
What_b = fn_What(mu_b)
if D / What_b > rhomax
error("Invalid state encountered at @ r = $(r)!")
end
if D < rhomin
# @info "Invalid state encountered at @ r = $(r)!"
# @info "Atmosphering ..."
# println((; D, rhomax, rhomin))
rho = rhoatm
ϵ = eos.cold_eos(InternalEnergy, Density, rho)
vr = 0.0
p = eos.cold_eos(Pressure, Density, rho)
# p = eos(Pressure, Density, rho, InternalEnergy, ϵ)
else
# do we have to adjust the bracket?
debug = false
# if D / What_b < rhomax < D
if rhomax < D
pblm_rhomax = Roots.ZeroProblem(mu -> D / fn_What(mu) - rhomax, (mu_a, mu_b))
mu_b = Roots.solve(pblm_rhomax, Roots.A42())
if isnan(mu_b)
error("Failed to correct upper interval bound")
end
debug = true
end
# if D / What_b < rhomin < D
if D / What_b < rhomin
mu_rng = collect(LinRange(mu_a, mu_b, 1000))
fn = @. D / fn_What(mu_rng) - rhomin
data = hcat(mu_rng, fn)
# display(data)
# display(D / What_b)
# display(D)
# display(rhomin)
pblm_rhomin = Roots.ZeroProblem(mu -> D / fn_What(mu) - rhomin, (mu_a, mu_b))
mu_b = Roots.solve(pblm_rhomin, Roots.A42())
if isnan(mu_b)
error("Failed to correct lower interval bound")
end
println("New interval: ($mu_a, $mu_b) @ r = $r")
debug = false
end
function master_fn(mu)
vhat = min(mu * s, v0)
What = 1 / sqrt(1 - vhat^2)
rhohat = D / What
ϵ_plus_one_over_What = What * (q - mu * s2)
ϵhat = ϵ_plus_one_over_What + vhat^2 * What^2 / (1 + What)
# enforce EOS bounds
_ϵmin, _ϵmax = fn_ϵmin(rhohat), fn_ϵmax(rhohat)
if rhohat < rhomin || ϵhat < _ϵmin
# TODO Should we reset What here?
rhohat, ϵhat = rhomin, _ϵmin
elseif rhohat > rhomax || ϵhat > _ϵmax
rhohat, ϵhat = rhomax, _ϵmax
end
# compute thermodynamic vars
phat = eos(Pressure, Density, rhohat, InternalEnergy, ϵhat)
ahat = phat / (rhohat * (1 + ϵhat))
hhat = eos(Enthalpy, Density, rhohat, InternalEnergy, ϵhat)
# compute master function
# extracted from (35), see hint in paragraph after (46)-(48)
nu_A = (1 + ahat) * ϵ_plus_one_over_What
nu_B = (1 + ahat) * (1 + q - mu * s2)
nu = max(nu_A, nu_B)
f = mu - 1 / (nu + s2 * mu)
return f
end
if debug
mu_rng = collect(LinRange(mu_a, mu_b, 1000))
fn = @. master_fn(mu_rng)
data = hcat(mu_rng, fn)
display(data)
end
mf_a = master_fn(mu_a)
mf_b = master_fn(mu_b)
if mf_a * mf_b >= 0
println((; D, Sr, τ, p, r, B, mf_a, mf_b, mu_a, mu_b))
error("Master function is ill conditioned, @ r = $(r)!")
end
pblm = Roots.ZeroProblem(master_fn, (mu_a, mu_b))
mu = Roots.solve(pblm, Roots.A42() #= =^= algorithm 748 =#)
if isnan(mu)
error("failed to find root")
end
# evaluate primitives
v = min(mu * s, v0)
W = 1 / sqrt(1 - v^2)
rho = D / W
ϵ_plus_one_over_W = W * (q - mu * s2)
ϵ = ϵ_plus_one_over_W + v^2 * W^2 / (1 + W)
# enforce EOS bounds
ϵmin, ϵmax = fn_ϵmin(rho), fn_ϵmax(rho)
if rho < rhomin || ϵ < ϵmin
rho, ϵ = rhomin, ϵmin
elseif rho > rhomax || ϵ > ϵmax
rho, ϵ = rhomax, ϵmax
end
# compute thermodynamic vars
p = eos(Pressure, Density, rho, InternalEnergy, ϵ)
vr = Sr*mu/D
end
@returns rho, vr, ϵ, p
end
#######################################################################
# Conservative variable fixing #
#######################################################################
function conservative_fixing(D, Sr, τ, B, equation)
@unpack eos, atm_density, atm_threshold = equation
rhomin = atm_density * atm_threshold
rhoatm = atm_density
if D < rhomin
D = rhoatm
end
if τ < 0
τ = 1e-12
end
# (48)-(50) in Deppe paper where we put B=0
# that = τ / D
# muhat = 0
# Solution to (52) when B=0
What = 1 + τ / D
# in spherical symmetry: S^2 = S_r^2 / B^2
# Hence, Smax^2/S^2 = Srmax^2 / Sr^2
S2 = Sr^2 / B^2
factor = sqrt( (1-1e-12) * (What^2-1) * D^2 / (S2 + D^2 * 1e-16) )
if factor > 1
Sr /= factor * (1 + 1e-8)
end
return D, Sr, τ
end