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.