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.

GRHD: Implement FV+slope limiters for doublecartoon version (!174)

Implements the FV+slope limiter machinery for the doublecartoon formulation.

While testing I noticed that the boundary conditions for a reflection boundary were implemented incorrectly.
To this end I quickly drafted a `BdryBufferedData` type with which one can abstract the buffer zones away.
This solution is not optimal performance wise, because it incurs quite some allocations (although I think it should be stack-allocateable...).

~~I have not bothered to update the spherical1d formulation, but instead inserted a `TODO`. So this needs to be updated in a separate PR.~~
Nevermind, I already fixed the todos, in order to make the reftests work.

Fixing the bdry conditions helped to improve the quality of the results, however, we are not even close yet to results from the literature.
Because of that I also now implemented the Roe method instead of the LLF method to compute the numerical fluxes, because that method
was also used in the paper Font et al 2000 PRD 61 044011.
This improves things too, but we are still worse than those results.
In particular, we use more resources and work in lower dimensions than that paper and yet we have larger errors ...
So this requires further investigations.

---

Also adds the following new parameters:
- `fv_numerical_flux`: the double cartoon version now implements the LLF and Roe flux.
- new option for `slope_limiter_method`: `MC`, this limiter is taken from Font et al 2002, PRD 65 084024.

SRHD: Fix FV slope limiter method to use up-to-date primitives in evolution (!173)

Due to the `cons2prim` methods being tied to `D, S, tau` registered as `dynamic_variables`, we need some juggling to get the primitives of the current step.

This inconvenience is due to
1. not (yet) separating the time stepping from the flux computation for the FV methods,
2. and the 'hard coding' of the primitives to the `cons2prim` methods.

Side note: We face the same issue with the `GRHD` version.

---

This PR also adds some more reftests to cover the `blastwave1` exercise.
To this end we also updated the FV solver to work with non-periodic bdry conditions with bdry values taken from initial data.

GRHD: Random improvements to TOV star test (!172)

- Implement AV treatment for `doublecartoon` formulation
- Implement Reconstruction treatment for `doublecartoon` formulation
- clean up logic in `callbacks.jl`
- Implement `fv_central` for `doublecartoon` formulation.

Also adds new parameters to `GRHD` section:
- `c2p_enforce_atm`: If true then atmosphere is enforced, regardless if its causal or not.
- `freeze_vars`: A list of variables for which `rhs` should be set to zero. Atm only supported by the doublecartoon formulation.

New parameter for the `TCI` section:
- `threshold_inv`: Invert the threshold check.

GRHD: Add initial data for Bondi accretion with infalling matter (!169)

This adds the following parameters to the `GRHD` section
- `id`: new option `bondi_accretion_infall`
- `id_bondi_accretion_infall_rs`: Radius below which the low-density fluid is initially set.
- `id_bondi_accretion_infall_ρl`: The density of the low-density fluid.
- `analyze_error_reference_solution`: the name of a reference solution that is to be used for the error estimation

---

This test is taken from Bugner's thesis 2017, and a paper by Radice, Rezolla 2011.

---

This PR also disables the questionable rescaling of the viscosity that was added for the TOV tests.
If needed it can be restored, but should be hidden behind a proper parameter.

CI: Add another CI target that runs the integration tests with 8 threads (!170)

@with_signature: Add `batch` option which enables multithreading of for loops using `Polyester.jl` (!168)

Usage:
```julia
@with_signature [batch=true] myfunction(eq::SomeEquation)
   @accepts ...
   # compute goes here
   @returns ...
end
```

I did some minimal profiling. There is a noticeable speed-up, but I did not measure how effective it is.
For the time being, I am happy when it just runs a little faster.

Evolve: use inline if in ode_step! for LSERK4 (!167)

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.