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/HRSC/mda.jl

@@ -41,7 +41,7 @@ function compute_viscosity!(
   # h = L / (K * N)
   h = L / K
 
-  u = reshape(u, layout(mesh))
+  u = dg1d.vreshape(u, layout(mesh))
   uhat = invV * u  # modal coefficients
 
   # perfect decay modal coefficients, added to enhance sensitivity to
@@ -65,19 +65,33 @@ function compute_viscosity!(
   rng_N = 1:N
   data_x = @. - log(rng_N)
 
+  mod_uh          = zeros(Float64, N)
+  skylined_mod_uh = zeros(Float64, N)
+  data_y          = zeros(Float64, N)
+
   for j in 1:K
     # (3.11)
     # only consider coefficients 2:Np (which corresponds to 1:N), because
     # the 1st (0th) one corresponds to a constant mode that is irrelevant
     # for the smoothness esimation.
-    mod_uh = @. sqrt( uhat[2:end,j]^2 + cellnorm2_u[j] * vhat^2 )
+    @. mod_uh = sqrt( uhat[2:end,j]^2 + cellnorm2_u[j] * vhat^2 )
     # (3.12)
     # max(1,N-1) inserted to avoid a 0:end range which can happen when only using two points
-    skylined_mod_uh = [ maximum(mod_uh[min(n,max(1,N-1)):end]) for n in 1:N ]
+    # skylined_mod_uh = [ maximum(mod_uh[min(n,max(1,N-1)):end]) for n in 1:N ]
+    for n = 1:N
+      skylined_mod_uh[n] = @views maximum(mod_uh[min(n,max(1,N-1)):end])
+    end
 
     # data points in log scale
-    data_y = @. log(skylined_mod_uh)
-    data_y[isinf.(data_y)] .= -20.0
+    # data_y = @. log(skylined_mod_uh)
+    @. data_y = log(skylined_mod_uh)
+    for n = 1:N
+      y = data_y[n]
+      if isinf(y)
+        data_y[n] = -20.0
+      end
+    end
+    # data_y[isinf.(data_y)] .= -20.0
     _, _, tau = dg1d.least_sqr_linear(data_x, data_y)
     isnan(tau) && error("undesired nan encountered")