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
@@ -300,5 +310,6 @@ end
 
   end
 
+  # @returns D, S, tau, rho, v, eps, p
   @returns rho, v, eps, p
 end