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/initialdata.jl

@@ -6,6 +6,9 @@ function initialdata!(env, P::Project, prms, mesh::Mesh1d)
   initialdata_hrsc!(env, P, P.hrsc, prms, mesh)
   initialdata_tci!(env, P, P.tci, prms, mesh)
   initialdata_bdrycond!(env, P, dg1d.DirichletBC2(), mesh)
+  if get(P.prms, :id_smooth, false)
+    initialdata_smooth!(env, P, mesh)
+  end
   return
 end
 
@@ -163,19 +166,16 @@ end
 function initialdata_equation_of_state(
     ::Val{:sod_shocktube}, eos::EquationOfState.IdealGas, mesh)
 
-  TODO()
-
+  # from EFL paper: arxiv 2202.08839
   @unpack x = get_static_variables(mesh.cache)
   (xmin, xmax), = mesh.extends
 
   @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
-  @toggled_assert eos.gamma == 5/3
+  @toggled_assert eos.gamma ≈ 1.4
   x0 = 0.0
-  # ρ0 = [ (xi < x0) ? 1.0  :  0.125  for xi in x ]
-  # v0 = [ (xi < x0) ? 0.0  :  0.0    for xi in x ]
-  # p0 = [ (xi < x0) ? 1.0  :  0.1    for xi in x ]
-
-  xl, xr = -0.5, 0.5
+  ρ0 = [ (xi < x0) ? 1.0  :  0.125  for xi in x ]
+  v0 = [ (xi < x0) ? 0.0  :  0.0    for xi in x ]
+  p0 = [ (xi < x0) ? 1.0  :  0.1    for xi in x ]
 
   (ρ0, v0, p0)
 end
@@ -225,7 +225,7 @@ function initialdata_equation_of_state(
   @unpack x = get_static_variables(mesh.cache)
   (xmin, xmax), = mesh.extends
 
-  @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+  @toggled_assert isapprox(xmin, -0.5) && isapprox(xmax, 0.5)
   @toggled_assert eos.gamma == 5/3
   x0 = 0.0
   ρ0 = [ (xi < x0) ? 10.0  : 1.0 for xi in x ]
@@ -242,7 +242,7 @@ function initialdata_equation_of_state(
   @unpack x = get_static_variables(mesh.cache)
   (xmin, xmax), = mesh.extends
 
-  @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+  @toggled_assert isapprox(xmin, -0.5) && isapprox(xmax, 0.5)
   @toggled_assert eos.gamma == 5/3
   x0 = 0.0
   ρ0 = [ (xi < x0) ? 1.0  : 1.0    for xi in x ]
@@ -430,7 +430,7 @@ function initialdata_hrsc!(env, P::Project, hrsc::Union{HRSC.EntropyViscosity,HR
 
   @unpack mu, cellmax_v = get_cell_variables(env.cache)
   @unpack max_v         = get_static_variables(env.cache)
-  @unpack D             = get_dynamic_variables(env.cache)
+  @unpack D, S, tau     = get_dynamic_variables(env.cache)
   dg1d.new_broadcast_volume!(maxspeed, P.equation, env.mesh)
   Npts, K = layout(env.mesh)
   mat_max_v = reshape(view(max_v, :), (Npts, K))
@@ -438,8 +438,8 @@ function initialdata_hrsc!(env, P::Project, hrsc::Union{HRSC.EntropyViscosity,HR
     cellmax_v[k] = maximum(view(mat_max_v, :, k))
   end
   HRSC.compute_viscosity!(
-      mu,
-      D, cellmax_v,
+      mu, deepcopy(mu),
+      (D,S,tau), cellmax_v,
       mda_hrsc.av)
 end
 
@@ -516,3 +516,26 @@ function initialdata_bdrycond!(env, P::Project2d, bdrycond::dg1d.DirichletBC2, m
   init_flx_tau_y .= flx_tau_y
   init_max_v     .= max_v
 end
+
+
+#######################################################################
+#                       Initial data smoothing                        #
+#######################################################################
+
+
+initialdata_smooth!(env, P::Project, mesh::Mesh{FVElement}) = nothing
+function initialdata_smooth!(env, P::Project, mesh)
+  !isregistered(mesh.cache, :mu) && return
+  @unpack D, S, tau = get_dynamic_variables(mesh.cache)
+  @unpack mu = get_cell_variables(mesh.cache)
+  bernstein = BernsteinReconstruction(mesh)
+  flag = zeros(mesh.K)
+  for k = 1:mesh.K
+    flag[k] = Float64(mu[k] > 0.0)
+  end
+  HRSC.reconstruct!(D,    flag, bernstein, isperiodic=true)
+  HRSC.reconstruct!(S,    flag, bernstein, isperiodic=true)
+  HRSC.reconstruct!(tau,  flag, bernstein, isperiodic=true)
+  dg1d.new_broadcast_volume!(cons2prim_kastaun, P.equation, mesh)
+  dg1d.new_broadcast_volume!(maxspeed, P.equation, mesh)
+end