remove comments

src/GRHD/cons2prim.jl

@@ -110,14 +110,10 @@ 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.ρ
-  # atm_threshold = equation.atmosphere.threshold
   ρmin = atm_density * atm_threshold
   @unpack K, Gamma = eos
   pmin = K * ρmin^Gamma
@@ -126,9 +122,7 @@ end
   ρatm = atm_density
   h0 = 1e-3
 
-  ####
-  # cons2prim from Kastaun et al. arXiv:2005.01821
-  ####
+  #### cons2prim from Kastaun et al. arXiv:2005.01821
 
   γuurr = 1/γrr
   q = τ / D # (22)-(23)
@@ -186,43 +180,27 @@ end
 
     # root bracket adjustment, cf. Appendix
     if D / What_b < ρmax < D
-      # @warn "adjusting upper root bracket"
-
       # 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
-
-      # 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"
-      # @show D/What_b, ρmin, D
-
       # 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
-
-      # @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)")
     end
 
     # infamous closure bug ...
@@ -259,20 +237,12 @@ end
     mf_a = master_fn(mu_a)
     mf_b = master_fn(mu_b)
     mf_ab = mf_a*mf_b
-    # 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)
       error("Failed to find root")
     end
 
-    #####
-    # evaluate primitives
-    #####
+    #### evaluate primitives
 
     # eq (36) [Kastaun] : h*W = 1/mu
     # eq (15) [Kastaun] : ρ*W = D
@@ -297,11 +267,8 @@ end
     v  = min(v, v0)
     v2 = v^2
     vr = sign(Sr) * sqrt(v2/γuurr)
-    # vr = sign(Sr) * min(abs(mu * rr), abs(v0))
 
     # eq (68) or (40)
-    # v = min(mu * r, v0)
-    # @show v2, vr, v0
     W = 1.0 / sqrt(1.0 - v2)
     # eq (41)
     ρ = D / W
@@ -312,12 +279,6 @@ end
 
     # enforce EOS bounds
 
-    # Deppe at all also have a atmopshere treatment for this
-    # if ρatm <= ρ <= ρmin
-    #   @assert v^2 < 1e-4
-    # end
-
-    # @assert ρ > ρmin ρ
     if ρ < ρmin
       # atmosphere II
 
@@ -342,37 +303,11 @@ end
       # error()
       # adjusting τ such that energy density is within bounds
       # cf. 3rd paragraph in sec. III.A in Kastaun paper
-      # @show ϵ, ϵmin, rcoord
       ϵ = ϵ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
 
-
-    # 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
-    #   # atmosphering
-    #   ρ, ϵ, v, W = max(ρmin, D), ϵmin, 0.0, 1.0
-    #   D = ρ
-    #   τ =
-    #   # ρ, ϵ, v, W = ρmin, ϵmin, 0.0, 1.0
-    # elseif ρ > ρmax || ϵ > ϵmax
-    #   @warn "ρmax || ϵmax reached"
-    #   flag_cons2prim_atm2 = 1
-    #   ρ, ϵ, v = ρmax, ϵmax, Si*mu/D
-    #   D = ρ
-    # end
-
     # compute thermodynamic vars
     p = eos(Pressure, Density, ρ, InternalEnergy, ϵ)
     if p < 0

src/GRHD/initialdata.jl

@@ -58,12 +58,6 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spheric
   equation.atm_factor = prms["atm_factor"]
   equation.atm_threshold = prms["atm_threshold_factor"]
 
-  # deps, ideps = newsignature(grhd_sources)
-  # deps  = renamevars(deps, Dict( (1,2,3,4) .=> (:r,:θ,:ϕ,:t) ))
-  # ideps = renamevars(ideps, Dict( (1,2,3,4) .=> (:r,:θ,:ϕ,:t) ))
-  # C = NamedTuple( (var,zeros(length(mesh))) for var in vcat(deps,ideps) )
-  # display(keys(C))
-
   # to avoid problems with 1/sin(θ) etc., this is also done in initialdata/SphericalAccretion.jl
   θ = π/2
   @unpack M = data