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)

GRHD: Revert changes from !114 regarding spherical1d formulation (!127, !127)

Although we had a reftest added for the TOV test with `spherical1d`, we had not enabled the comparison of data. As such the test was almost useless. This is fixed here for this test. But we should add more reftests for other formulations in separate PRs.

GRHD Fixup atmosphere freezing logic (!126, !126)

We need to account for the fact that global solution is degenerated on cell interfaces.
Because of that we need to cast our data stencil wider when we check for resets on a bdry.

GRHD: Implement reflection bdry conditions for doublecartoon formulation and TOV test. (!125, !125)

I decided to implement the symmetry conditions now differently than I did for the other formulations (`spherical1d, rescaled_spherical1d, valencia1d`).
That is, I now compute the bdry values of the state variable separately before I call `bdryllf_...`.
The values are computed such that
- for variables that are even wrt `r -> -r` trafos, the x derivative of the nodal approximation at the origin vanishes,
- for variables that are odd wrt `r -> -r` trafos, the value of the nodal approximation at the origin vanishes.
Given this data (stored in `bdry_D, bdry_Sx, bdry_tau`), I then impose the bdry conditions weakly by using the LLF flux and the previously computed
bdry values as exterior states.
I also tested computing the exterior states according to `q^+ = -q^- + 2 q_BC`, where `q_BC` are the previously computed values, but did not find that
this has significant effect.

What is interesting is that this implementation (doublecartoon + reflection bdry conds.) can now be evolved without symmetrizing the state vector explicitly
and it appears to be stable, e.g. no blowup at the origin.
I wonder whether adapting this way of imposing bdry conditions for other formulations would also help with stability, although I am not sure anymore
if they suffered from blowups in the reflected case.
Note that all other formulations suffer from the problem of division by zero, which we dynamically fix up by explicit symmetrization. Imposing the bdry
conditions then differently won't fix that issue.

However, for now we keep the bdry conditions in the other formulations as they currently are.

---

More thoughts about bdry values:
- Like I said above, all other formulations use the values from initial data as bdry values, which makes sense because the solution is supposed to be static.
One should check if using instead bdry values such that the symmetry conditions are satisfied actually works better and/or is factually more accurate.
The latter might be plausible, because numerical dissipation (espcially with AV) might just cause a decrease in the central density.
Using instead the initial data as bdry values would try to keep the central density as it is, despite the star solution being altered.
- For the mirrored domain test we can have a grid setup in which the cells faces fall onto the origin. It would make sense to test if using similarly
constructed interface values for the LLF flux at the interfaces that lie at the origin. Maybe that could also cure the stability issue we see in
all the mirrored runs.

* update parameter docstrings

* add minimal description for each formulation

* remove dead code

* split equation.jl into spherical1d.jl, rescaled_spherical1d.jl, valencia1d.jl, doublecartoon.jl

* manually undo the insertion of the refvol factor into spherical1d
formulation

the formulae without the refvol factor should already be correct
in particular, they correspond to the choice where all three
metrics (physical one, conformal one and reference one) all agree
I will document this properly in a follow up PR were I also
disentangle all formulations

* rename llf and bdryllf to refer to spherical1d formulation

* update example parfile

* fix derivatives of induced metric in initial data

* remove unused variables

* fix sign in 4 metric determinant

* fixup minus and fix doublecartoon corrections at x=0

* fix todo conflict

* some todos

* implement symmetrizer around origin, but disable it for now

* fix computation of state variables in llf_doublecartoon

* some debug statements

* implement the doublecartoon correction

* implement the doublecartoon correction

* remove not evolved y,z components and fixup nomenclature

* add doublecartoon implementation

* factor out data reader from initial data setup

* update parfile

* add more tov examples

* fix spherical1d to account for reference volume element

* fixup impose_symmetry! for valencia1d

* reset tau when it falls below 0

this is needed to run with lower atmosphere densities

* impose symmetry on state variables before c2p

* callbacks: fix atm freezing logic

* c2p: todos

* implement valancia1d version that works with tov test

* add grhd_bondi example that uses valencia1d formulation

* implement valencia1d version

* add valencia1d option for formulation parameter

* c2p: more verbose log output

* register more variables for rescaled_spherical1d

* fix dispatch in setup.jl

* verify formulation parameter

* add bdryllf for rescaled_spherical1d

* stop gap for unhandled cases of singularities in impose_symmetry functions

* add simulation time to progress meter output

* implement cartoon derivative expressions for mixed rank (1,1) tensor

* add docstrings

* fix imports

* add cartoon macros

* update 2d reftests to check 2d output

also reduced number of output variables and
times in order to save storage

* IntegrationTests: error when reference data is missing

* fix filename of 2d output

* also test 0d output with fv_muscl_burgers_sine reftest

* fix u_norm2 computation for FVElement

* reftests: update advection_sine_analysis

* implement u_norm2 0d analysis variable