GRHD: Implement reflection bdry conditions for doublecartoon formulation and TOV test.

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.