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

@@ -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