adjust tau when energy density exceeds lower bound

src/SRHD/cons2prim.jl

@@ -98,7 +98,7 @@ end
   rhomax = Inf
   rhoatm = atm_density
   # epsmin adjusted for IdealGas
-  epsmin, epsmax = 1e-10, Inf
+  epsmin, epsmax = 0.0, Inf
   h0 = 1e-3
 
   B = 1.0 # in special relativity
@@ -264,16 +264,16 @@ end
     # regularize velocity according to Dumbser et al arxiv:2307.06629
     # eqs (82)-(84)
     y, y0 = rhohW2, 1e-4
-    if y > y0
+    # if y > y0
       # v = S / mu
       # compute v accoriding to eq (68) or (40) [Kastaun]
       v = sign(Si) * min(mu * r, 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
+    # 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
 
     # eq (68) or (40)
     # v = min(mu * r, v0)
@@ -292,24 +292,39 @@ end
     #   @assert v^2 < 1e-4
     # end
 
-    if rho < rhomin || eps < epsmin
-      # @warn "atmosphering II @ x = $x ($rho, $eps)"
-      # atmosphering
-      rho, eps, v, W = max(rhomin, D), epsmin, 0.0, 1.0
-      # rho, eps, v, W = rhomin, epsmin, 0.0, 1.0
-    elseif rho > rhomax || eps > epsmax
-      @warn "rhomax || epsmax reached"
-      rho, eps, v = rhomax, epsmax, Si*mu/D
+    @assert rho > rhomin
+    if eps < epsmin
+      # adjusting tau such that energy density is within bounds
+      # cf. 3rd paragraph in sec. III.A in Kastaun paper
+      eps = epsmin
+      taub4 = tau
+      # invert this for tau: eps = W * (q - mu * r2) + v^2 * W^2 / (1.0 + W), q = tau / D
+      tau = (eps / W  - v^2 * W / (1.0 + W) + mu * r2) * D
+      @warn "atmosphering II @ x = $x : (taub4,tau) = ($taub4, $tau)"
     end
+    # if rho < rhomin || eps < epsmin
+    #   @warn "atmosphering II @ x = $x : (rho,eps) = ($rho, $eps)"
+    #   flag_cons2prim_atm2 = 1
+    #   # atmosphering
+    #   rho, eps, v, W = max(rhomin, D), epsmin, 0.0, 1.0
+    #   D = rho
+    #   tau =
+    #   # rho, eps, v, W = rhomin, epsmin, 0.0, 1.0
+    # elseif rho > rhomax || eps > epsmax
+    #   @warn "rhomax || epsmax reached"
+    #   flag_cons2prim_atm2 = 1
+    #   rho, eps, v = rhomax, epsmax, Si*mu/D
+    #   D = rho
+    # end
 
     # compute thermodynamic vars
     p = eos(Pressure, Density, rho, InternalEnergy, eps)
-    if p <= 0
+    if p < 0
       error((rho,eps))
     end
 
   end
 
   # @returns D, S, tau, rho, v, eps, p
-  @returns rho, v, eps, p
+  @returns tau, rho, v, eps, p
 end