ScalarEq: Implement modal Bspline2 DG evolution (!213)

Similar to !207, but using the new 2nd order Bspline basis.

Currently supports only bare evolutions and no HRSC etc.

Also: This is a quasi-nodal Ansatz similar to the Bernstein Ansatz, i.e. it is not interpolating, but the modal coefficients are related to some sampling of the data on the grid.
This means that one has to sum the modal coefficients onto the grid to obtain the nodal values, which can be done with `interpolate` or a Vandermonde matrix.
Atm this is not consistently done in the ScalarEq project, hence, care needs to be examined when interpreting results!

This adds multiple things:
- A `Bspline` module for working with 2nd order Bsplines.
- Extends `DGElement` with the `kind=:modal_bspline2` option.
- Updates docstrings and refactors some tests.
- Adds a `LG` module for Legendre-Gauss quadrature.
- Adds a plotting script `plot/bspline.jl` to tinker with the `Bspline2` approximation.

Also adds a new option to the `Mesh` parameters section:
- `kind = "modal_bspline2"`: Utilize a 2nd order modal Bspline Ansatz for the DG evolution.

CI: try to reuse precompilation cache at least for IntegrationTests (!212)

dg1d: Fix convention for Bernstein polynomials (!211)

We now establish the convention that all basis functions are indexed by one.

Also had to update the `modal_bernstein_advection_sine` reftest data, because !201 skipped running CI on main and me cancelling some other CI runs.

We require that all PRs are rebased before merging, so running them twice does not help anything.
This should free up the CI queue.

ScalarEq: Implement modal Bernstein DG evolution (!207)

MWE to check consistency of operators.

This is not intended to work with any HRSC schemes or 2d etc.

dg1d: Rename `SpectralElement` to `DGElement` and generalize its purpose (!206)

Generalization should allow it represent DG operators that are not necessarily based on nodal approximations via Lagrange polynomials.

Atm we implement the `nodal_bernstein` version, which is actually pseudo-nodal and not yet tested. Will be done in a follow-up PR.

Along the way we also implement some operators for Bernstein and Legendre bases.

tests: add kw only_failed to `runtests` command (!205)

Unify interface for LGL and GLGL rules and update `SpectralElement` accordingly:
- most routine should just work with one argument `Npts`
- add more doc strings
- remove unused methods `linear_dependency_matrix` and `expansion_matrix`
- remove unused fields of `SpectralElement`

GRHD: fix another edge case where c2p_set_atmosphere_on_failure should overrule (!203)

dg1d: substep callbacks have to be evolved after updating the solution (!202)

dg1d: update `lgl` and `glgl` integration methods (!201)

- put the integrand as the first argument
- add integration methods that can work with a function argument
- add doc strings
- update usage of `integrate` functions in rest of code

GRHD: fix atmosphere freezing and evolution logic (!199)

Fixes the following issues with the current implementation:
- `update_atm_domain_of_dependence!` did not distinguish between points
  1. whose neighbors are also in the bulk,
  2. where one neighbor is on an interface,
  3. where the point itself is on an interface.

  These cases are already properly handled in 2d.
- `stop_atmosphere_evolution_*` methods did not account for the `c2p_freeze_atm` flag, which essentially prevented the atmosphere to ever evolve.
- `update_atm_domain_of_dependence!` was activated by setting `c2p_dynamical_atm`, but it should be controlled by `c2p_enforce_causal_atm`



---

Additional fixes with less impact:
- remove unused parameter `id_atmosphere_density` from `GRHD` section
- `c2p_set_atmosphere_on_failure` should really suppress any errors appearing in the cons2prim

I am surprised that I did not need this before.

This is intentionally a hacky implementation as I just need that for the thesis right now.

---

This also implements
- the total rest-mass computation for the GRHD project as a first application,
- adds helper functions for the `Output` struct to read all data at once.

Fixup to !193

GRHD: add tov symmetrization for all DG rhs' that use the spherical1d (!195)

dg1d: make `main` accept paths to output directories (!196)

dg1d: do not strip trailing newline chars as this might corrupt git patches (!194)

dg1d:  allow parameter overriding on cmd line (!193)

Allows to use `main(parfile, "Mesh.n" => "1")` to override the value given in the parfile.

---

To do this we
- parse the TOML file
- flatten the dict-of-dicts to a single dict with the keys joined by a `.`
- insert the parameters
- unflatten the dict to a dict-of-dicts
- write this to the output directory

Unfortunately, `TOML.jl` does not allow to parse into `OrderedDict`.
This causes the newly written parfiles to be messed up.
Furthermore, by manually re-writing the parfiles we also do not preserve any comments or formatting.

dg1d: Avoid generating nested output folders in case parameter file comes from inside an output folder (!192)

GRHD: add toggle to allow AV to sense based on the log10 of the variable (!190)

Adds a new parameter to the `GRHD` section:
- `av_sensor_logabs_D`: If `true` apply the AV sensor to `log(|D|)`, otherwise use `D`, where `D` is the state variable for the rest-mass density (including the volume factor).

Evolve: align stepper logic between all ERK implementations (!189)

This assumes that the first row in the `a` matrix of the Butcher tableau for ERK methods has only zeros, which seems to be the case for the all the methods have atm.

---

Problem explanation:

So far the LSERK4 and the remaining ERK steppers worked slightly different.
The former used for the first substep the initial data array `u0`, and `u!`  for the remaining steps. That shortcut avoids a copy and appears natural for the LSERK4 algorithm.
OTOH, in the ERK stepper we first copy `u0` into `u!`, then accumulate some substep `k`s into `u!`, and then call the RHS with `u!`. (the reason we always copied `u0` over was that we did not account for the above assumption on `a` before)

The problem was now that in the GRHD evolutions we sometimes alter the state vector in the RHS computation.
This influenced the results to the extent that we could not run the (newly added) reftest only with the LSERK4 method.
What was happening was that when using the ERK method we never altered `u0` directly, only `u!`, whereas `u0` was altered in case the GRHD RHS decided to do so on the first step.
Subsequent substeps of the ERK method depend upon `u0`, and so this broke things.

GRHD: ban `ParameterDB` usage from `rhs` due to type instabilities (!188)

GRHD: enable `bc = from_id` for dynamic Bondi accretion test (!187)

ScalarEq: adds the av_recompute_substeps option (!186)

In analogy to how its done for SRHD.

Adds the following parameter to the `ScalarEq` section:
- `av_recompute_substeps`: If `true` then all callbacks related to `HRSC` are ran in substeps of the time stepper.

ScalarEq: Fix time step computation for DG method (!185)

similar !183

ScalarEq: Implement general slope limiter methods (!184)

Reuses the implementation from `SRHD`.

The following parameters are added/updated to the `ScalarEq` section:
- `hrsc` (updated): added option `slope_limiter`
- `slope_limiter_method` (new): similar to `SRHD.slope_limiter_method`
- `slope_limiter_method_tvb_M` (new): similar to `SRHD.slope_limiter_method_tvb_M`

SRHD: Use correct timestep scaling for vanilla DG (!183)

Evolution: fix substep callbacks shall be called after rhs evaluation (!182)

I haven't used substep callbacks much, so this should not have an effect on current experiments.

GRHD: Make AV method work with sperical1d (!181)

Move some methods to were they belong.
And also enable AV method to work with the spherical1d formulation.

GRHD: Add option to disable the cartoon derivatives when running with the cartoon formulation (!180)

Adds the following option to the `GRHD` section:
- `doublecartoon_disable_derivatives`: If `true`, then the doublecartoon derivative terms are not added to the source terms.

It appears that this option is senseless, because why even use the doublecartoon formulation in the first place.
However, it turns out that these derivative terms severely degrade the quality of any TOV star tests I tried so far,
and removing those terms from the equation helps a ton. The reason is likely that these terms contain 1/r (or 1/x) factors
which will amplify any noise coming from the TOV star's surface and propagating inwards.

Note that removing those terms changes the PDEs we solve.
In fact, removing these terms was already considered in Bugner's paper Physical Review D, 2016.

GRHD: more specific boundry conditions [only incorporated for FV+slope_limiter] (!179)

This adds new options under `GRHD.bc`:
- `tov_reflective_origin`: impose reflection symmetry at coordinate origin, impose vacuum values at outer bdry
- `tov_symmetric_domain`: impose vacuum values at both domain boundary ends

Also streamlined the setup of the bdry values for FV.
While doing so I realized that the inflow bdry in the Bondi test was wrongly treated.
There we cannot extrapolate interior data, nor just copy what we have in the interior, unless we are very far out (which in my test we aren't).
Instead, what we do now there is use the initial data, which is the analytical solution and extrapolate that.
This improves the behavior of the solution to no longer drift after the first time step. Nice.

The formulae isn't complicated, but I still messed up doing so.

Fixing it makes TOV evolutions in symmetric domains now remain symmetric, kind of.

Because `minmod` is used with almost all slope limiters, I had to update quite a few reftest too, although they overall quality of the results did not change much.

I wondering whether this fix makes now the WENO version of DG also work, will check in a separate PR.

GRHD: fix implementation of `atm_evolve` (!177)

Follow up to !176

Adds the following new parameter:
- `GRHD.atm_evolve`: If `false` then the RHSs of all state variables are set to zero at any grid point where we identified the data to describe atmosphere.

This is parameter is implemented to work with all formulations and schemes, hopefully.

GHRD: Add Marquina's numerical flux for FV+slope limiter and doublecartoon (!175)

Adds the following new parameters
- `GRHD.fv_numerical_flux = "marquina"`: Numerical flux implementation following Font et al 2000, PRD 61 044011; see Font et al 2002, PRD 65 084024 for TOV tests.

We also update the TOV initial data module to include a function with which one should be able to generate all the configurations that we use in the ref tests.