cons2prim: use non-zero values for rho, eps thresholds

src/GRHD/cons2prim.jl

@@ -119,10 +119,11 @@ end
   # atm_density = equation.atmosphere.ρ
   # atm_threshold = equation.atmosphere.threshold
   ρmin = atm_density * atm_threshold
+  @unpack K, Gamma = eos
+  pmin = K * ρmin^Gamma
+  ϵmin = pmin / ((Gamma-1)*ρmin)
   ρmax = Inf
   ρatm = atm_density
-  # ϵmin adjusted for IdealGas
-  ϵmin, ϵmax = 0.0, Inf
   h0 = 1e-3
 
   ####
@@ -341,6 +342,7 @@ 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