Integrated and discretised right-hand side (SSInfall)

This directory hosts the original methodology as discussed in

@article{10.1093/mnras/stac1991,
    author = {Zhang, Han and Weinzierl, Tobias and Schulz, Holger and Li, Baojiu},
    title = "{Spherical accretion of collisional gas in modified gravity I: self-similar solutions and a new cosmological hydrodynamical code}",
    journal = {Monthly Notices of the Royal Astronomical Society},
    volume = {515},
    number = {2},
    pages = {2464-2482},
    year = {2022},
    month = {07},
    abstract = "{The spherical collapse scenario has great importance in cosmology since it captures several crucial aspects of structure formation. The presence of self-similar solutions in the Einstein-de Sitter (EdS) model greatly simplifies its analysis, making it a powerful tool to gain valuable insights into the real and more complicated physical processes involved in galaxy formation. While there has been a large body of research to incorporate various additional physical processes into spherical collapse, the effect of modified gravity (MG) models, which are popular alternatives to the Λ cold dark matter paradigm to explain the cosmic acceleration, is still not well understood in this scenario. In this paper, we study the spherical accretion of collisional gas in a particular MG model, which is a rare case that also admits self-similar solutions. The model displays interesting behaviours caused by the enhanced gravity and a screening mechanism. Despite the strong effects of MG, we find that its self-similar solution agrees well with that of the EdS model. These results are used to assess a new cosmological hydrodynamical code for spherical collapse simulations introduced here, which is based on the hyperbolic partial differential equation engine ExaHyPE 2. Its good agreement with the theoretical predictions confirms the reliability of this code in modelling astrophysical processes in spherical collapse. We will use this code to study the evolution of gas in more realistic MG models in future work.}",
    issn = {0035-8711},
    doi = {10.1093/mnras/stac1991},
    url = {https://doi.org/10.1093/mnras/stac1991},
    eprint = {https://academic.oup.com/mnras/article-pdf/515/2/2464/45192162/stac1991.pdf},
}

A preprint of this document is available from https://arxiv.org/abs/2203.01174.

The present algorithm is rather simple. For the modified right-hand side, we need a solution to the Poisson equation \( - \Delta \phi = -\pi G/\rho \). Actually, it is not the integral over \( \rho \) but over the difference of \( \rho \) to the baseline density. We have a fundamental solution to this operator and hence can compute \( \phi \) (or \( \nabla \phi \) respectively, which is the quantity of interest), but for this, we to know \( \int -\pi G/\rho dx \).

We discretise time and make the right-hand side of the current time step feed into the solution of the subsequent time step. Our algorithm reads as follows:

1. Evolve the Euler equations subject to a known right-hand side.
2. Compute \( \int -\pi G/\rho dx \) of the new solution.
3. Prepare a new right-hand side for the next time step.

The integral \( \int -\pi G/\rho dx \) can be complex. However, we know that the function has rotation symmetry. Therefore, we discretise the domain into layers similar to an onion. We bucket the solution. In the example above, we have four layers \( \rho _1, \rho _2, \rho _3, \rho _4 \). That is, we write \( \int -\pi G/\rho dx = \sum _{i} \int _{r_i,r_i+1} -\pi G/\rho V_i dr \), where r is the radius away from the centre, and V is the volume of the corresponding shell. Every time we have updated a cell, we look into which onion ring the cell centre falls and accumulate the corresponding \( \rho _i \). We discretise the onion, and then we sort cells into these shells. This gives us an approximation for the volume, where the term \( V_i \) is implicitly approximated through the cells' volumes.

Realisation

The logic of the accumulation of the overdensities is realised within the class applications::exahype2::euler::sphericalaccretion::MassAccumulator. The accumulator is not only able to keep track of the discretised density distribution within shells, it also provides the tools to interpolate between those buckets.

The outcome feeds into the Euler equation's source term via a gravitional model such as applications::exahype2::euler::sphericalaccretion::addGravitationalSource_AlphaCDM(). It accepts the mass resulting from the overdensity in a region and translates it into a vector driving the PDE.

Workflow

The integration of the modified source term into the PDE is straightforward, as ExaHyPE's API offers a source term callback function already.

To realise the accumulation, we add a volume- or dof-wise postprocessing kernel to the solution. All ExaHyPE solvers offer such a feature (the exahype2.solvers.rkdg.actionsets.DoFWisePostprocessSolution action set is an example).

Finally, we have to ensure that the mass accumulator's MPI exchange routines are called after each and every time step.