utilize c2p_threshold_velocity_limiter parameter in cons2prim

src/GRHD/cons2prim.jl

@@ -135,7 +135,7 @@ end
   verbose = false
   c2p_reset_ϵ, c2p_init_admissible = 0.0, 1.0
   # setup some constants
-  @unpack atm_density, atm_threshold, eos = equation
+  @unpack atm_density, atm_threshold, c2p_threshold_velocity_limiter, eos = equation
   ρmin = atm_density * atm_threshold
   @unpack K, Gamma = eos
   pmin = K * ρmin^Gamma
@@ -143,7 +143,8 @@ end
   ρmax = Inf
   ρatm = atm_density
   h0 = 1e-3
-  y0 = 1e-4
+  y0 = c2p_threshold_velocity_limiter
+  e0 = y0^2/2 # fudge factor for filtered division in velocity limiter
   Dmin = ρmin * 1 # W = 1
   # assuming diagonal metric
   γuuxx = 1/γxx
@@ -215,7 +216,7 @@ end
   #   @warn "S2max check failed at $xcoord, $zcoord because S2 = $S2 vs S2max = $S2max"
   # end
 
-  limit_v = let D=D, v0=v0, Sx=Sx, Sz=Sz, r=r, γuuxx=γuuxx, γuuzz=γuuzz, _y0=y0
+  limit_v = let D=D, v0=v0, Sx=Sx, Sz=Sz, r=r, γuuxx=γuuxx, γuuzz=γuuzz, _y0=y0, _e0=e0
     # infamous closure bug ...
     # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
     mu -> begin
@@ -226,7 +227,7 @@ end
         min(mu*r, v0)
       else
         _fy = 2.0 * (_y/_y0)^3 - 3.0 * (_y/_y0)^2 + 1.0
-        _f = _y / (_y^2 + _fy * 5e-9)
+        _f = _y / (_y^2 + _fy * _e0)
         # no need to use rx*D here, because D,rx,Sx should remain unchanged throughout
         _vx = Sx * _f
         _vz = Sz * _f
@@ -257,7 +258,7 @@ end
   end
 
   if (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && D < ρmin)
-    @warn D ρmin xcoord zcoord c2p_freeze_atm c2p_reset_atm
+    @warn "This should not have happened" D Sx Sz τ ρmin xcoord zcoord c2p_freeze_atm c2p_reset_atm
     TODO("Upsi")
   elseif (c2p_freeze_atm == 0.0 && D < ρmin) || (c2p_freeze_atm > 0.0 && c2p_reset_atm > 0.0)
     @assert eos isa EquationOfState.Polytrope
@@ -359,7 +360,7 @@ end
     else
       c2p_limit_vr = 1.0
       fy = 2.0 * (y/y0)^3 - 3.0 * (y/y0)^2 + 1.0
-      f = y / (y^2 + fy * 5e-9)
+      f = y / (y^2 + fy * e0)
       vx = Sx * f
       vz = Sz * f
     end