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

@@ -109,7 +109,7 @@ end
   r = abs(ri)
   v0 = abs(v0i)
 
-  if v0 >= 1
+  if v0 > 1
     error((v0,z0i,Si,D))
   end
 
@@ -134,10 +134,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