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

@@ -6,9 +6,9 @@ timestep(env, P, hrsc) = timestep(env, P, hrsc, env.mesh)
 
 
 function compute_maxspeed(mesh, equation, cache)
-  @unpack v = get_static_variables(cache)
+  @unpack max_v = get_static_variables(cache)
   dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
-  vmax = dg1d.absolute_maximum(v)
+  vmax = dg1d.absolute_maximum(max_v)
   vmax_limit = 1e4
   if vmax > vmax_limit
     @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
@@ -31,15 +31,15 @@ dtfactor(mesh::Mesh1d{SpectralElement}) = mesh.element.N
 function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, mesh::Mesh1d)
   @unpack cache, mesh = env
   @unpack equation    = P
+  @unpack max_v = get_static_variables(cache)
 
   # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
-  dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
-  dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
-  max_v = get_variable(cache, :max_v)
+  # dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+  # dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
   vmax = dg1d.absolute_maximum(view(max_v,:))
-  vmax_limit = 0.99
+  vmax_limit = 0.9999999
   if vmax > vmax_limit
-    @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
+    @warn "Limiting timestep due to maximum speed $vmax exceeding $vmax_limit"
     vmax = vmax_limit
   end
   smoothed_mu = get_variable(cache, :smoothed_mu)
@@ -206,11 +206,53 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction,
 end
 
 
+function fd_diff_1d!(_df, _f, _x, mesh)
+  @unpack Npts = mesh.element
+  for (df, f, x) in zip(eachcell(mesh,_df),eachcell(mesh,_f),eachcell(mesh,_x))
+    for n in 2:Npts-1
+      fl,fm,fr = f[n-1],f[n],f[n+1]
+      xl,xm,xr = x[n-1],x[n],x[n+1]
+      dll = (xm-xr)/((xl-xm)*(xl-xr))
+      dlm = (2*xm-xl-xr)/((xm-xl)*(xm-xr))
+      dlr = (xm-xl)/((xr-xl)*(xr-xm))
+      df[n] = fl*dll + fm*dlm + fr*dlr
+    end
+    fl,fm,fr = f[1],f[2],f[3]
+    xl,xm,xr = x[1],x[2],x[3]
+    dll = (2*xl-xr-xm)/((xl-xm)*(xl-xr))
+    dlm = (xl-xr)/((xm-xl)*(xm-xr))
+    dlr = (xl-xm)/((xr-xl)*(xr-xm))
+    df[1] = fl*dll + fm*dlm + fr*dlr
+    fl,fm,fr = f[Npts-2],f[Npts-1],f[Npts]
+    xl,xm,xr = x[Npts-2],x[Npts-1],x[Npts]
+    dll = (xr-xm)/((xl-xm)*(xl-xr))
+    dlm = (xr-xl)/((xm-xl)*(xm-xr))
+    dlr = (2*xr-xl-xm)/((xr-xl)*(xr-xm))
+    df[Npts] = fl*dll + fm*dlm + fr*dlr
+  end
+end
+
+
 function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
     bdryconds, ldg_bdryconds, av_bdryconds, ::Mesh1d)
 
-  @unpack cache, mesh           = env
-  @unpack equation              = P
+  @unpack cache, mesh = env
+  @unpack equation, prms = P
+  reg = prms.av_regularization
+
+  if reg === :mono
+    rhs_mono!(cache, mesh, equation, P)
+  elseif reg === :covariant
+    rhs_covariant!(cache, mesh, equation, P)
+  else
+    TODO(reg)
+  end
+
+end
+
+
+function rhs_mono!(cache, mesh, equation, P)
+
   @unpack D, S, tau             = get_dynamic_variables(cache)
   @unpack flx_D, flx_S, flx_tau,
           ldg_D, ldg_S, ldg_tau,
@@ -224,6 +266,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
           bdry_ldg_D, bdry_ldg_S,
           bdry_ldg_tau,
           bdry_smoothed_mu      = get_bdry_variables(cache)
+  @unpack x = get_static_variables(cache)
 
   ## solve auxiliary equation: q + ∂x u = 0
   dg1d.new_broadcast_faces!(ldg_nflux, equation, mesh)
@@ -232,7 +275,12 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
   compute_rhs_weak_form!(ldg_D,   D,   nflx_D,   mesh)
   compute_rhs_weak_form!(ldg_S,   S,   nflx_S,   mesh)
   compute_rhs_weak_form!(ldg_tau, tau, nflx_tau, mesh)
+  # fd_diff_1d!(ldg_D,   D,   x, mesh)
+  # fd_diff_1d!(ldg_S,   S,   x, mesh)
+  # fd_diff_1d!(ldg_tau, tau, x, mesh)
+
 
+  # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
   dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
   dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
 
@@ -270,6 +318,85 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
 end
 
 
+function rhs_covariant!(cache, mesh, equation, P)
+
+  @unpack D, S, tau             = get_dynamic_variables(cache)
+  @unpack flx_D, flx_S, flx_tau,
+          ldg_D, ldg_S, ldg_tau,
+          smoothed_mu, max_v,
+          rho, v, eps, p        = get_static_variables(cache)
+  @unpack rhs_D, rhs_S, rhs_tau = get_rhs_variables(cache)
+  @unpack nflx_D, nflx_S, nflx_tau,
+          bdry_D, bdry_S, bdry_tau,
+          bdry_rho, bdry_v, bdry_eps,
+          bdry_p, bdry_max_v,
+          bdry_ldg_D, bdry_ldg_S, bdry_ldg_tau,
+          bdry_rhs_D, bdry_rhs_S, bdry_rhs_tau,
+          bdry_smoothed_mu      = get_bdry_variables(cache)
+  @unpack x = get_static_variables(cache)
+
+  # compute rhs for ∂t u + ∂x f(x) = 0
+
+  # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
+  dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+  dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, D,     bdry_D)
+  dg1d.interpolate_face_data!(mesh, S,     bdry_S)
+  dg1d.interpolate_face_data!(mesh, tau,   bdry_tau)
+  dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
+  dg1d.interpolate_face_data!(mesh, rho,   bdry_rho)
+  dg1d.interpolate_face_data!(mesh, v,     bdry_v)
+  dg1d.interpolate_face_data!(mesh, eps,   bdry_eps)
+  dg1d.interpolate_face_data!(mesh, p,     bdry_p)
+
+  dg1d.new_broadcast_volume!(flux, equation, mesh)
+
+  # numerical fluxes
+  dg1d.new_broadcast_faces!(llf,    equation, mesh)
+  dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
+
+  compute_rhs_weak_form!(rhs_D,   flx_D,   nflx_D,   mesh)
+  compute_rhs_weak_form!(rhs_S,   flx_S,   nflx_S,   mesh)
+  compute_rhs_weak_form!(rhs_tau, flx_tau, nflx_tau, mesh)
+
+  # now compute rhs for ∂t u + ∂x f(x) = ∂x( μ (-∂t+∂x) u )
+
+  ## solve auxiliary equation: q + ∂x u = 0
+  dg1d.new_broadcast_faces!(ldg_nflux,  equation, mesh)
+  dg1d.new_broadcast_bdry!(ldg_bdryllf, equation, mesh)
+  compute_rhs_weak_form!(ldg_D,   D,   nflx_D,   mesh)
+  compute_rhs_weak_form!(ldg_S,   S,   nflx_S,   mesh)
+  compute_rhs_weak_form!(ldg_tau, tau, nflx_tau, mesh)
+  # fd_diff_1d!(ldg_D,   D,   x, mesh)
+  # fd_diff_1d!(ldg_S,   S,   x, mesh)
+  # fd_diff_1d!(ldg_tau, tau, x, mesh)
+
+  # need to recompute, because ldg computation overwrites nflx_{D,S,tau}
+  # numerical fluxes
+  dg1d.new_broadcast_faces!(llf,    equation, mesh)
+  dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
+
+  dg1d.new_broadcast_volume!(av_flux_covariant, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, ldg_D,       bdry_ldg_D)
+  dg1d.interpolate_face_data!(mesh, ldg_S,       bdry_ldg_S)
+  dg1d.interpolate_face_data!(mesh, ldg_tau,     bdry_ldg_tau)
+  dg1d.interpolate_face_data!(mesh, rhs_D,       bdry_rhs_D)
+  dg1d.interpolate_face_data!(mesh, rhs_S,       bdry_rhs_S)
+  dg1d.interpolate_face_data!(mesh, rhs_tau,     bdry_rhs_tau)
+  dg1d.interpolate_face_data!(mesh, smoothed_mu, bdry_smoothed_mu)
+
+  dg1d.new_broadcast_faces!(av_nflux_covariant, equation, mesh)
+
+  compute_rhs_weak_form!(rhs_D,   flx_D,   nflx_D,   mesh)
+  compute_rhs_weak_form!(rhs_S,   flx_S,   nflx_S,   mesh)
+  compute_rhs_weak_form!(rhs_tau, flx_tau, nflx_tau, mesh)
+
+  return
+end
+
+
 #######################################################################
 #                                 2d                                  #
 #######################################################################