make shock tube and blast wave tests run for SRHD when using MDA AV

The key step was to adapt a 'covariant' regularization that looks like

\partial_t D + \partial_x f_D^x = - \partial_x ( \mu ( \partial_t - \partial_x ) D ),

where the time derivative on the RHS is new compared to previous attempts. We approximate \partial_t D by the - \partial_x f_D^x.

Note that this is not truely covariant, but (\partial_t - \partial_x) ( \mu ( \partial_t - \partial_x ) D) would be. Problem is that we don't know how to approximate the second time derivatives or the mixed time-space derivatives. But it works this way, so leave it for now.

This PR adds example parameter files that allow to reproduce the SRHD tests from https://arxiv.org/pdf/2202.08839.pdf. To run the blastwave2 test it was important to add correct the evolved variable tau whenever the energy density eps turned slightly negative, as suggested in https://arxiv.org/pdf/2005.01821.pdf III.A. What is interesting is that atm we don't need an atmosphere description, but that's probably because those tests do not examine low densities.

---

Side note: Realized that !75 (merged) promoted Mesh to a parameteric type and now most of the HRSC and TCI structs suffer from type instabilities. Need to fix this in a separate PR.

src/SRHD/cons2prim.jl

@@ -11,7 +11,7 @@
   atm_threshold = equation.atmosphere.threshold
   rhomin = atm_density * atm_threshold
 
-  B = 1.0 # in special relativity
+  # B = 1.0 # in special relativity
 
   #####
   # conservative variable fixing following Deppe et al. arXiv:2109.12033
@@ -19,17 +19,23 @@
 
   fixed = false
 
-  if D < rhomin
-    # @warn "fixing D $D @ x = $x"
+  if D < rhomin || tau < 0.0
     D = rhomin
-    fixed = true
-  end
-
-  if tau < 0.0
-    # @warn "fixing tau $tau @ x = $x"
     tau = 1e-12
     fixed = true
   end
+
+  # if D < 0.001 * rhomin
+  #   # @warn "fixing D $D @ x = $x"
+  #   D = rhomin
+  #   fixed = true
+  # end
+  #
+  # if tau < 0.0
+  #   # @warn "fixing tau $tau @ x = $x"
+  #   tau = 1e-12
+  #   fixed = true
+  # end
   # fix tau such that v < 0.999
   # if tau < 0.0
   #   @warn "fixing tau $tau @ x = $x"
@@ -40,33 +46,37 @@
   #   fixed = true
   # end
 
-  # (48)-(50) in arXiv:2109.12033 where we put B=0
-  # that = tau / D
-  # muhat = 0
-  # Solution to (52) when B=0
-  # What = 1.0 + tau / D
+  # (48)-(50) in arXiv:2109.12033 we put B=0, ̂μ = 0
+  # solution to (52) is then
+  What = 1.0 + tau / D
+
+  S2 = S^2
+  D2 = D^2
+  ϵs = 1e-12
+  factor = sqrt( (1.0-ϵs) * (What^2-1.0) * D2 / (S2 + D2 * 1e-16) )
+  # @show factor, What
+  if factor < 1.0 && What > 1.0
+    # @warn "fixing ..."
+    Sold = S
+    S = 0.99 * sign(S) * sqrt( (1.0-ϵs) * (What^2-1.0) * D2 - D2 * 1e-16 )
+    # @warn "S  = $Sold => $S @ x = $x"
+    S2 = S^2
+    factor = max(1.0, sqrt( (1.0-ϵs) * (What^2-1.0) * D2 / (S2 + D2 * 1e-16) ))
+    # @show factor
+    @assert factor >= 1.0
+    fixed = true
+  end
 
+  # h0 = 1e-3
   # S2 = S^2 / B^2
-  # factor = sqrt( (1.0-1e-12) * (What^2-1.0) * D^2 / (S2 + D^2 * 1e-16) )
-  # if factor <= 1.0
-  #   @warn "fixing S $S @ x = $x"
-  #   S = sign(S) * sqrt( (1.0-1e-12) * (What^2-1.0) * D^2 - D^2 * 5e-12  )
-  #   S2 = S^2 / B^2
-  #   factor = sqrt( (1.0-1e-12) * (What^2-1.0) * D^2 / (S2) )
-  #   @assert factor >= 1.0
-  #   fixed = true
+  # r2 = S2 / D^2
+  # z02 = r2 / h0^2
+  # v02 = z02 / (1.0 + z02)
+  # max_v2 = 0.999
+  # if v02 > max_v2
+  #   S = sign(S) * sqrt(max_v2/(1-max_v2) * h0^2 * D^2)
   # end
 
-  h0 = 1e-3
-  S2 = S^2 / B^2
-  r2 = S2 / D^2
-  z02 = r2 / h0^2
-  v02 = z02 / (1.0 + z02)
-  max_v2 = 0.999
-  if v02 > max_v2
-    S = sign(S) * sqrt(max_v2/(1-max_v2) * h0^2 * D^2)
-  end
-
   flag_consfixed = Float64(fixed)
 
   @returns D, S, tau, flag_consfixed
@@ -78,7 +88,7 @@ end
 #######################################################################
 
 
-@with_signature [legacy=false] function cons2prim_kastaun(equation::AbstractEquation)
+@with_signature [simd=false,legacy=false] function cons2prim_kastaun(equation::AbstractEquation)
   @accepts D, S, tau, x
 
   @unpack eos = equation
@@ -88,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
@@ -109,7 +119,7 @@ end
   r = abs(ri)
   v0 = abs(v0i)
 
-  if v0 >= 1
+  if v0 > 1
     error((v0,z0i,Si,D))
   end
 
@@ -134,10 +144,31 @@ end
   if D < rhomin
     @warn "atmosphering I"
 
+    error((D,S,tau))
+
+    @assert eos isa EquationOfState.IdealGas
+    @unpack Gamma = eos
+
+    # employ atmosphere values by taking the zero temperature limit
+    # for an ideal gas eos, we have that eos.Gamma becomes the Gamma for
+    # a polytrope, see grhd notes for derivation
+
     rho = rhoatm
-    eps = eos.cold_eos(InternalEnergy, Density, rho)
     v = 0.0
-    p = eos.cold_eos(Pressure, Density, rho)
+    W = 1.0
+    # p = K * rho^Gamma
+    # eps = Gamma * rho^
+    # D = W * rho
+    # S = rho * (1+eps+p/rho) * W^2 * v
+    # tau = rho * (1+eps+p/rho) * W^2 - p - D
+
+    # if eos isa ColdEquationOfState
+    #   eps = eos.cold_eos(InternalEnergy, Density, rho)
+    #   p = eos.cold_eos(Pressure, Density, rho)
+    # else
+    #   eps = eos(InternalEnergy, Density, rho)
+    #   p = eos(Pressure, Density, rho)
+    # end
 
   else
 
@@ -233,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)
@@ -261,23 +292,38 @@ 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
     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 rho, v, eps, p
+  # @returns D, S, tau, rho, v, eps, p
+  @returns tau, rho, v, eps, p
 end