19.3 Self-gravity

The self-consistent gravity algorithm supplied with FLASH computes the Newtonian gravitational field produced by the matter. The produced potential function satisfies Poisson's equation (19.2). This unit's implementation can also return the acceleration field ${\bf g}({\bf x})$ computed by finite-differencing the potential using the expressions

$$g_{x;ijk} = {1\over2\Delta x}\left(\phi_{i-1,j,k} - \phi_{i+1,j,k}\right) + {\cal O}(\Delta x^2)$$

$$g_{y;ijk} = {1\over2\Delta y}\left(\phi_{i,j-1,k} - \phi_{i,j+1,k}\right) + {\cal O}(\Delta y^2)$$ (19.3)

$$g_{z;ijk} = {1\over2\Delta z}\left(\phi_{i,j,k-1} - \phi_{i,j,k+1}\right) + {\cal O}(\Delta z^2) .$$

In order to preserve the second-order accuracy of these expressions at jumps in grid refinement, it is important to use quadratic interpolants when filling guard cells at such locations. Otherwise, the truncation error of the interpolants will produce unphysical forces at these block boundaries.

Two algorithms are available for solving the Poisson equations: Gravity/GravityMain/Multipole and Gravity/GravityMain/Multigrid. The initialization routines for these algorithms are contained in the Gravity unit, but the actual implementations are contained below the Grid unit due to code architecture constraints.

The multipole-based solver described in Sec:GridSolversMultipole for self gravity is appropriate for spherical or nearly-spherical mass distributions with isolated boundary conditions. For non-spherical mass distributions higher order moments of the solver must be used. Note that storage and CPU costs scale roughly as the square of number of moments used, so it is best to use this solver only for nearly spherical matter distributions.

The multigrid solver described in Sec:GridSolversMultigrid is appropriate for general mass distributions and can solve problems with more general boundary conditions.

The tree solver described in 8.10.2.4 is appropriate for general mass distributions and can solve problems with both isolated and periodic boundary conditions set independently in individual directions.

19.3.1 Coupling Gravity with Hydrodynamics

The gravitational field couples to the Euler equations only through the momentum and energy equations. If we define the total energy density as

$$\rho E \equiv {1\over 2}\rho v^2 + \rho\epsilon ,$$ (19.4)

where $\epsilon$ is the specific internal energy, then the gravitational source terms for the momentum and energy equations are $\rho{\bf g}$ and $\rho{\bf v}\cdot{\bf g}$, respectively. Because of the variety of ways in which different hydrodynamics schemes treat these source terms, the gravity module only supplies the potential $\phi$ and acceleration ${\bf g}$, leaving the implementation of the fluid coupling to the hydrodynamics module. Finite-difference and finite-volume hydrodynamic schemes apply the source terms in their advection steps, sometimes at multiple intermediate timesteps and sometimes using staggered meshes for vector quantities like ${\bf v}$ and ${\bf g}$.

For example, the PPM algorithm supplied with FLASH uses the following update steps to obtain the momentum and energy in cell $i$ at timestep $n+1$

$$(\rho v)_i^{n+1} = (\rho v)_i^n + {\Delta t\over 2} g_i^{n+1} \left(\rho_i^n + \rho_i^{n+1}\right)$$

$$(\rho E)_i^{n+1} = (\rho E)_i^n + {\Delta t\over 4} g_i^{n+1} \left(\rho_i^n + \rho_i^{n+1}\right)\left(v_i^n + v_i^{n+1}\right) .$$ (19.5)

Here $g_i^{n+1}$ is obtained by extrapolation from $\phi_i^{n-1}$ and $\phi_i^n$. The ${\tt Poisson}$ gravity implementation supplies a mesh variable to contain the potential from the previous timestep; future releases of FLASH may permit the storage of several time levels of this quantity for hydrodynamics algorithms that require more steps. Currently, ${\bf g}$ is computed at cell centers.

Note that finite-volume schemes do not retain explicit conservation of momentum and energy when gravity source terms are added. Godunov schemes such as PPM, require an additional step in order to preserve second-order time accuracy. The gravitational acceleration component $g_i$ is fitted by interpolants along with the other state variables, and these interpolants are used to construct characteristic-averaged values of ${\bf g}$ in each cell. The velocity states $v_{L,i+1/2}$ and $v_{R,i+1/2}$, which are used as inputs to the Riemann problem solver, are then corrected to account for the acceleration using the following expressions

$$v_{L,i+1/2} \rightarrow v_{L,i+1/2} + {\Delta t\over 4}\left(g^+_{L,i+1/2} + g^-_{L,i+1/2}\right)$$

$$v_{R,i+1/2} \rightarrow v_{R,i+1/2} + {\Delta t\over 4}\left(g^+_{R,i+1/2} + g^-_{R,i+1/2}\right) .$$ (19.6)

Here $g^\pm_{X,i+1/2}$ is the acceleration averaged using the interpolant on the $X$ side of the interface ($X=L,R$) for $v\pm c$ characteristics, which bring material to the interface between cells $i$ and $i+1$ during the timestep.

19.3.2 Tree Gravity

The Tree implementation of the gravity unit in physics/Gravity/GravityMain/Poisson/BHTree is meant to be used together with the tree solver implementation Grid/GridSolvers/BHTree/Wunsch. It either calculates the gravitational potential field which is subsequently differentiated in subroutine Gravity_accelOneRow to obtain the gravitational acceleration, or it calculates the gravitational acceleration directly. The latter approach is more accurate, because the error due to numerical differentiation is avoided, however, it consumes more memory for storing three components of the gravitational acceleration. The direct acceleration calculation can be switched on by specifying bhtreeAcc=1 as a command line argument of the setup script.

The gravity unit provides subroutines for building and walking the tree called by the tree solver. In this version, only monopole moments (node masses) are used for the potential/acceleration calculation. It also defines new multipole acceptance criteria (MACs) that estimate the error in gravitational acceleration of a contribution of a single node to the potential (hereafter partial error) much better than purely geometrical MAC defined in the tree solver. They are: (1) the approximate partial error (APE), and (2) the maximum partial error (MPE). The first one is based on an assumption that the partial error is proportional to the multipole moment of the node. The node is accepted for calculation of the potential if

$$D^{m+2} > \frac{GMS_\mathrm{node}^m}{\Delta a_\mathrm{p,APE}}$$ (19.7)

where $D$ is distance between the point-of-calculation and the node, $M$ is the node mass, $S_\mathrm{node}$ is the node size, $m$ is a degree of the multipole approximation and $\Delta a_\mathrm{p,APE}$ is the requested maximum error in acceleration (controlled by runtime parameter grv_bhAccErr). Since only monopole moments are used for the potential calculation, the most reasonable choice of $m$ seems to be $m=2$. This MAC is similar to the one used in Gadget2 (see Springel, 2005, MNRAS, 364, 1105).

The second MAC (maximum partial error, MPE) calculates the error in acceleration of a single node contribution $\Delta a_\mathrm{p,MPE}$ according to formula 9 from Salmon&Warren (1994; see this paper for details):

$$\Delta a_\mathrm{p,MPE} \le \frac{1}{D^2}\frac{1}{(1-S_\mathrm{no... ...}\left( \frac{3\lceil B_2\rceil}{D^2} - \frac{2\lfloor B_3 \rfloor}{D^3}\right)$$ (19.8)

where $B_n = \Sigma_i m_i \vert\mathbf{r}_i-\mathbf{r}_0\vert^n$ where $m_i$ and $\mathbf{r}_i$ are masses and positions of individual grid cells within the node and $\mathbf{r}_0$ is the node mass center position. Moment $B_2$ can be easily determined during the tree build, moment $B_3$ can be estimated as $B_3^2 \ge B_2^3/M$. The maximum allowed partial error in gravitational acceleration is controlled by runtime parameters grv_bhAccErr and grv_bhUseRelAccErr (see 19.4.1).

During the tree walk, subroutine Gravity_bhNodeContrib adds contributions of tree nodes to the gravitational potential or acceleration fields. In case of the potential, the contribution is

$$\Phi = -\frac{GM}{\vert\vec{r}\vert}$$ (19.9)

if isolated boundary conditions are used, or

$$\Phi = -GMf_\mathrm{EF,\Phi}(\vec{r})$$ (19.10)

if periodic periodic or mixed boundary conditions are used. In case of the acceleration, the contributions are

$$\vec{a}_g = \frac{GM\vec{r}}{\vert\vec{r}\vert^3}$$ (19.11)

for isolated boundary conditions, or

$$\vec{a}_g = GMf_\mathrm{EF,a}(\vec{r})$$ (19.12)

for periodic or mixed boundary conditions. In In the above formulae, $G$ is the constant of gravity, $M$ is the node mass, $\vec{r}$ is the position vector between point-of-calculation and the node mass center and $f_\mathrm{EF,\Phi}$ and $f_\mathrm{EF,a}$ are the Ewald fields for the potential and the acceleration (see below).

Boundary conditions can be isolated or periodic, set independently for each direction. If they are periodic at least in one direction, the Ewald method is used for the potential calculation (Ewald, P. P., 1921, Ann. Phys. 64, 253). The original Ewald method is an efficient method for computing gravitational field for problems with periodic boundary conditions in three directions. Ewald speeded up evaluation of the gravitational potential by splitting it into two parts, $Gm/r=Gm \mathrm{erf}(\alpha r)/r + Gm \mathrm{erfc}(\alpha r)/r$ ($\alpha$ is an arbitrary constant) and then by applying Poisson summation formula on erfc terms, gravitational field at position $\vec{r}$ can be written in the form

$$\phi (\vec{r}) = -G \sum_{a=1}^N m_a \left( \sum_{i_1,i_2,i_3} A_... ...i_3}) \right) = -G \sum_{a=1}^N m_a f_\mathrm{EF,\Phi}(\vec{r}_a - \vec{r}) ,$$ (19.13)

the first sum runs over whole computational domain, where at position $\vec{r}_a$ is mass $m_a$. Second sum runs over all neigbouring computational domains, which are at positions $\vec{l}_{i_1,i_2,i_3}$ and $A_S(\vec{r},\vec{r}_a,\vec{l}_{i_1,i_2,i_3})$ and $A_L(\vec{r},\vec{r}_a,\vec{l}_{i_1,i_2,i_3})$ are short and long-range contributions, respectively. It is sufficient to take into account only few terms in eq. 19.13. The Ewald field for the acceleration, $f_\mathrm{EF,a}$, is obtained using a similar decomposition. We modified Ewald method for problems with periodic boundary conditions in two directions and isolated boundary conditions in the third direction and for problems with periodic boundary conditions in one direction and isolated in two directions.

The gravity unit allows also to use a static external gravitational field read from file. In this version, the field can be either spherically symmetric or planar being only a function of the z-coordinate. The external field file is a text file containing three columns of numbers representing the coordinate, the potential and the acceleration. The coordinate is the radial distance or z-distance from the center of the external field given by runtime parameters. The external field if mapped to a grid using a linear interpolation each time the gravitational acceleration is calculated (in subroutine Gravity_accelOneRow).