this helps the profile outputs to associate the measurements to the same line

- Don't impose initial data values at interior boundary. This is an inflow boundary, but we cheat around that and just use the interior flux there.
- Add a higher resolved version of the Bondi accretion data. This is necessary to get acceptable results in a convergence study.

SRHD: perform c2p when wrapping data in callback (!164)

TBH this is a bandaid to ensure that primitives are up-to-date before outputting.

In the future there should be a proper mechanism to ensure that.

GRHD: Use WeakBoundaries() for CubicSplineInterpolation (!163)

Sometimes our grid bounds don't align exactly with the bounds from the data. This boundary treatment is more forgiving in that sense.

Adds two new parameters:
- `variables0d_analysis_error`: A list of variable names for which the error wrt the initial data values and computed as the L2 norm shall be computed. For each enqueued variable name a `Output.variables0d` entry is added which consists of the variable name and a `_err_norm` postfix. E.g. `rho` adds a 0d output `rho_err_norm`.
- `variables1d_analysis_error`: Similar to `variables0d_analysis_error`, but instead of computing the L2 norm, the difference over the mesh is provided as an `Output.variables1d` entry with name equal to variable name plus `_err` postfix.

---

This is a rushed implementation, because I need that output now.

But this mechanism should be 1) partly factored out into the `dg1d` project, and 2) we should not need extra parameters.
About 2): It should be enough to ask for `rho_err` or `rho_err_norm` in 1d or 0d output fields. Doing so would require us to reserve names ending with `_err` and `_err_norm`, but that should be fine.

Furthermore, we should rename the postfixes to something like `_id_err` and `_id_err_l2`, to make their nature more clear.

EulerEq: Allow use of more than one AV indicator (!161)

This adds two new options
- `av_indicator_variables`: a list of variables to which the AV indicator
  should be applied to. This generates one viscosity per variable.
- `av_combine_operator`: an operater with which multiple viscosities
  are reduced to a single one.

We also re-activated the `euler_sod_shock_tube_entropyav` reftest, which now seems to give acceptable results.
Also updated the `euler2d_kelvin_helmholtz_entropyav` reftest, although that should not have been necessary.
However, 2d experiments are not declared stable anyways, so will need to check those later again regardless.

EulerEq: implement the `bc="from_id"` option for slope limited FV methods (!160)

EulerEq: use the correct initial data for the Sod shock tube test (!159)

GRHD: fix limiting of velocity components (!158)

This was a subtle mistake where applying the velocity limiter build into the Kastaun scheme was not correctly applied to components of the velocity.

EulerEq: Rework Navier-Stokes regularization (!155)

The goal is to reproduce plots for the shock tube tests from the Zingan+ 2013 paper.

We rework the Navier-Stokes regularization to work more like in the Zingan+ 2013 paper. Before we tried to mimic a Guermond+ 2014 paper, but they did not provide any shock tube tests.

This PR adds new parameters to the `EulerEq` section:
- `av_navierstoke_prho, av_navierstoke_pt`: these are effective Prandtl numbers that can be used to control how much of the viscosity coefficient for the momentum regularization should be added to the density and energy regularization.
- `id_smooth`: If true smooth initial data profile with a Bernstein filter
- `av_drag`: Dragging factor to throttle AV evolution
- `av_recompute_substeps`: If true AV computation is scheduled to run on each substep.

We also fixed a sever issue in one of the formulae.
In the computation of `maxspeed` we accidentally compute `q/p` instead of `q/rho` as one part of the wave speed.
Fixing this seems to have a positive influence.
I am curious how this will effect the ref tests.
Answer: I had to update all associated reftests, but all changes were minor.

Initial data configuration taken from Zingan+ 2013 and Toro 1999.

c2p: Add option `c2p_set_atmosphere_on_failure` (!153)

Adds a new option `c2p_set_atmosphere_on_failure` to the `GRHD` and `SRHD` sections.
If `true` and the c2p solver fails to invert the master function, we reset all variables to atmosphere.

We also fixed that the running with `verbose` in `c2p` no longer resest to atmosphere.

fix invocation of substep callbacks for explicit schemes (!151)

- Mainly refactors `c2p` in order to allow to ignore the total energy density in the recovery procedure.
This is possible for cold equations of state, where we don't need to evolve tau.

- The `c2p` now also incorporates a conservative fixing that we adapted from from arXiv:1306.4953, Appendix.

- We had to update reftest data for the SRHD simple wave and GRHD bondi tests.
Especially the GRHD bondi tests improved from not using tau in the reconstruction.

- We also updated the default value `atm_treshold_factor` so that it is now slightly large than one.
This fixes some errors where we misidentified atmosphere states as being below atmosphere, before the c2p procedure started.

- Because this PR was cherry-picked from the fa/doublecartoon+av branch, we also brought over the changes in the naming of variables of the doublecartoon formulation (which should have went into another PR, but I was too lazy to entangle those changes first).

----

I noticed that the number of allocs in `rhs` sections are now quite increased.
This is particularly bad for the SRHD tests. My guess is this is due to this refactor now passing on both EoSs as parameters to the c2p routine.

This PR replaces the SRHD's own c2p with that of GRHD.

The following parameters are added to the SRHD section:
- `id_smooth = false`: If `true` initial data is smoothed with a Bernstein filter.
- `av_regularization = "none"`: Controls which AV regularization is to be used.
- `c2p_cold_eos_parameters = [ 100.0, 2.0 ]`: Parameters of cold polytropic EoS for c2p, analogous to GRHD.
- `av_drag = 0.0`: AV dragging factor to control how rapid AV can change, analogous to GRHD and EulerEq.
- `av_recompute_substeps = false`: If `true` AV recomputation is triggered on every substep in the evolution.

Copy current Manifest.toml to output folder (!148)

GRHD: Use a dedicated cold EoS as a fallback in c2p (!147)

This adds new parameters to the `GRHD` section:
- `c2p_cold_eos_parameters`: Contains parameters for the `Polytrope` EoS that is used as a fallback in c2p.

This kind of 'abstraction' is needed so that we can re-purpose this c2p for the SRHD project.

EulerEq: Fix 2d equations and timestep scaling (!146)

- The 2d version did compute the flux incorrectly.
- We used too strong of a scaling for the time step, e.g. N^-3 instead of N^-2. This is fixed for 1d and 2d case.
- Updated `euler2d_isentropic_flow` ref test to no longer run into a shock, because that example doesn't use any HRSC.
- Updated `euler_sod_shock_tube` test to run with low resolution. Also shorten `tend` so that data does not reach boundary.

This addresses two issues:
- Utilize cold EoS when we limit velocity. This should avoid artificial spikes in epsilon (internal energy density) due to the limiting, because the derived relations used in the Kastaun paper no longer hold.
- Make `atm_threshold_factor = 1.0` the default. In some tests I noticed that having this greater than 1.0 an artificial plateau in the evolved density appears, which I did not like. This seems to fix that.

This adds and modifies the following parameters in the `GRHD` section:
- (new) `c2p_dyamic_atm`: Controls whether a fluid-atmosphere interface is allowed to evolve.
- (renamed) `c2p_freeze_atm_reset -> c2p_enforce_causal_atm`; also updated the doc string

With `c2p_dyamic_atm=true` we now recompute the `c2p_atm_reset` mask before determining the `c2p_freeze_atm` mask.
Previously, we only ever ran with the `false` option, which would mean that we artificially imposed the interface as it was given in the initial data.

GRHD: remove stale isadmissible check from cons2prim (!144)

It sems that atm we can get away with slightly negative tau values.
Furthermore, I think the logic for correcting tau wasn't really sound, and I could not retrace anymore how I obtained that fix.

GRHD: More fixes to the atmosphere freezing mask in 2d (!141)

- fixed two neighbor indices
- fixed next-neighbor check to also account for -1.0 (invalid) values

GRHD: Make parameters of cons2prim's velocity filter configurable (!140)

This adds the following parameter to the `GRHD` section:
- `c2p_threshold_velocity_limiter`: Threshold `y₀ = ρhW²` below which the velocity is limited by a filtered division following Dumbser+ 2024.

I found emperically that this value strongly influences the reflection and generation of perturbations in the TOV star tests.
In particular, using 1e-4 gives a PSD of the central density evolution in which all peak frequencies are wrong compared to literature (e.g. Three-dimensional numerical general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars, Font+, 2002, table 1).
Using instead a value of 1e-6 gives better results.
However, this seems to come at the expense of stability, because with lower limiter values we appear to more easily apply surface heating, e.g. epsilon or tau display artifacts at the surface.

The default value of `c2p_threshold_velocity_limiter` is now set to `1e-6` because of the above observation.
This also explains why updated the ref test in this PR.

We did not account for the fact that there a two
data points on an interface as well as two neighbors
(one in each direction).

This adds new parameters to the `Mesh` section:
- `filter_enable`: Enable spectral filtering to be applied after each full time step.
- `filter_parameters`: Parameters for exponential spectral filter, cf. (5.16) in Hesthaven, Warburton 2007.

GRHD: also impose EoS limits when computing final primitives (!137)

This went unnoticed so far, because we have been running only with polytropes and there epsilon is basically unused.

Evolve: Fix t argument for callback_substep for ode_step (!136)

Readd previously removed `load_parameters` (!135)

GRHD: Implement TOV test in 2d axisymmetry using the cartoon formulae (!134, !134)

The `examples/grhd_tov/dg_cartoon.toml` example can run for a good while without immediately crashing. This means that the bdry conditions (reflection symmetry across `x = 0` and `z = 0` axes) aren't so wrong and the surface treatment is doing something.

However, ... I realized that I have been running with smaller time steps than other schemes do. That is, mine scales as `~ N^{-3}`, but it should more be like `~ N^{-2}`.
This explains why I could run with comparatively large CFL numbers some of the 1d tests.
This should be investigated in a separate PR.
Changing to `~N^{-2}` here makes the cartoon test crash shortly after start, because of an issue in the cons2prim recovery. So perhaps we can't get away without any regularization scheme and I have to implement AV too.

Some fixes to dg1d.jl core (!133, !133)

- Adds `setup_env(parfile)` with which you can now directly construct `env, project, bdryconds, prms` from a parfile. These objects are the basic ingredients for an evolution.
- Fixes overspecialization in caching code.
- Fixes a bug where showing `CallbackTrigger` caused an error.

GRHD: Make `cons2prim` work with 2D velocities (!132, !132)

I am glad I have the reftests in place, because this was very delicate to implement correctly.
Part of this struggle was also to not introduce extra allocations from boxing of the closures, due to Julia sharing variables between scopes ...

This PR also applies a general clean up to the `cons2prim` function.

`GRHD`: Implement reflection bdry conditions for `valencia1d, spherical1d, rescaled_spherical1d` formulations (!131, !131)

It works nicely with `spherical1d, rescaled_spherical1d`, but fails miserably with `valencia1d`.
Not sure what the problem there is, because for both `rescaled_spherical1d` and `valencia1d` I need to interpolate onto the origin after unscaling the conservatives,
because the volume factor comes with an extra `r` that makes the conservatives vanish there.
However, for `valencia1d` and `spherical1d` I also have singular source terms, which I don't have in `rescaled_spherical1d`.

`@with_signature`: avoid tuple limit by using `get_variable` unpacking directly in function body (!130, !130)

The implementation of the `doublecartoon` formulation in `GRHD` finally hit the tuple length limit, above which tuples are stack allocated and, thus, incur a drastic performance hit due to extra allocs.

Here we update `@with_signature` so that we no longer have to pass arguments through tuples,
but instead we look them up with the `get_variable` method from `mesh.cache`.
This way we also do the unpacking now directly inside the macro-generated function body.

Its not an ideal solution, because now `@with_signature` can no longer be used without the `get_variable` method from `dg1d`, which increases coupling.
But for the moment I don't bother, and I shall worry about that whenever I decided to move `@with_signature` into a separate pkg, if ever.

---

The `get_variable` method and `Cache` could also use a bit of benchmarking and optimizations.
E.g. atm we perform three lookups
- `isregistered` to check if the variable is available
- `Cache.variables_to_fields` to determine in which variable group the variable can be found
- `getproperty(::Variable, ::Symbol)` to get the variable

Clearly, this could be improved, but we leave that for the future.

---

Obviously, I also updated the tests. But I refused to setup now `Cache`s for each test case.
Instead I just added a dispatch `get_variable(::NamedTuple, ::Symbol)` and a `cache` argument to the expanded function so that the tests didn't so much adjustments.

GRHD: Fix rescaled_spherical1d formulation (!128, !128)

A few issues have silently accumulated because we did not add a reftest when we introduced this formulation.
This PR fixes them and adds a reftest.

GRHD: Add reftest for doublecartoon formulation (!129, !129)