AN ADAPTIVE PARTICLE-MESH GRAVITY SOLVER FOR ENZO

The default gravity solver implemented in Enzo is based on a combination of Fourier and multigrid algorithms. We recall here the basics of this method, but for more details, see Bryan et al. (2014). First the gravitating potential is computed on the root grid in Fourier space using a fast Fourier transform (FFT, HE87). For periodic boundary conditions, the Green's function is generated directly in Fourier-space. For isolated boundary conditions, it is calculated first in real space to obtain the correct zero-padding, and transformed in the Fourier domain (James 1977). In both cases, the resulting potential is then transformed to the real domain and differentiated in order to obtain the face-centered (cell-centered) acceleration field depending on whether the hydro solver requires face-centered or cell-centered accelerations. On the subgrids, boundary conditions for the potential are obtained by interpolation from the parent grid and then Equation (1) is solved using a Gauss-Seidel multigrid relaxation method with Dirichlet boundary conditions (HE87). As we will highlight in Section 3, such algorithm may introduce errors on the refined grids due to inaccuracies in the coarse-grid solution. We should also emphasize that the particular multigrid solver described here is not strictly speaking a pure multigrid solver in the traditional sense. Indeed, classic multigrid solvers calculate the solution on all levels at the same time.

The adaptive particle-mesh solver is based on the particle-particle adaptive particle-mesh method (HE87, Couchman 1991). The general idea of the algorithm is to split the gravitational force between a long-range component, and one or more short-range components which are nonzero only for a narrow range of wavenumbers. The algorithm as described in this section is suitable for three-dimensional problems only, but could be extended in the future to treat one- and two-dimensional systems.

The calculation on the root grid is nearly identical to the one performed with the multigrid solver (although the Green's function is slightly modified). On a refined grid one first calculates the short-range component of the force. Therefore, the Green's function needs to be modified in order to account for the contribution from the smallest scales only. We thus need a smoothing function, which here is the sphere with uniformly decreasing density S(r):

$$\begin{equation} S(r) = \left\lbrace \begin{array}{@{} ll}\frac{48}{\pi a^4} \left(\frac{a}{2} - r \right) & {\rm if} \quad r \le a/2 \\ \\ 0 & {\rm otherwise,} \end{array} \right. \end{equation} \tag{ 3 }$$

where r is the radius and a is a positive parameter. The corresponding Fourier transform is

$$\begin{equation} \hat{S}_l (k) = \frac{12}{\eta ^4} \left(2 - 2\cos {\eta } - \eta \sin {\eta }\right), \end{equation} \tag{ 4 }$$

where k = |k| and η = ka/2. For each refinement level we take a = 3.4δ_l, where δ_l is the linear size of a grid cell on the current level l. The effects of the smoothing function and the smoothing parameter have been studied extensively in HE87, who concluded that this profile gives a better accuracy in three-dimensional schemes and this chosen value for a minimizes the numerical noise.

Minimizing the error in the mesh force leads to the optimal Green's function (HE87, Equation (8–22)):

$$\begin{equation} \hat{G}(\mathbf {k}) = \frac{\mathbf {\hat{D}} (\mathbf {k}) \cdot \sum _\mathbf {n} \hat{U}^2(\mathbf {k_n}) \mathbf {\hat{R}}(\mathbf {k_n})}{|\mathbf {\hat{D}} (\mathbf {k})|^2 \left[ \sum _\mathbf {n} \hat{U}^2(\mathbf {k_n}) \right]^2 }, \end{equation} \tag{ 5 }$$

where

$$\begin{equation} \mathbf {k_n} = \mathbf {k} +\mathbf {n} 2\pi /\delta _l, \quad \mathbf {n} \in [-\infty,+\infty ]^3 \end{equation} \tag{ 6 }$$

are the Brillouin zones, $\mathbf {\hat{D}}$ is the gradient operator, U is the assignment function, and R is the reference force. For the gradient operator we use either a direct Fourier-space solver ($\hat{D}(\mathbf {k}) = -i \mathbf {k}$) or a finite-difference representation (see HE87 for an explicit expression). For TSC deposition, the assignment function U is (HE87, Equations (8–42), (8–45))

$$\begin{equation} \hat{U}(\mathbf {k}) = \prod _{i=1}^3\, {\rm sinc}^3\left({\frac{k_i\delta _l}{2}} \right) \end{equation} \tag{ 7 }$$

and verifies

$$\begin{equation} \sum _\mathbf {n} \hat{U}^2(\mathbf {k_n}) = \prod _{i=1}^3 \left(1-\sin ^2{\frac{k_i\delta _l}{2}} + \frac{2}{15} \sin ^4{\frac{k_i\delta _l}{2}} \right). \end{equation} \tag{ 8 }$$

Finally, the reference force is linked to the smoothing function by

$$\begin{equation} \hat{\mathbf {R}}(\mathbf {k}) = -i \frac{\hat{S}_l^2 - \hat{S}_{l-1}^2 }{k^2} \mathbf {k}. \end{equation} \tag{ 9 }$$

$\hat{S_l}$ signifies the smoothing function with a cell size on the level $l\mathbf\; \ge\; 1$. For the coarsest level (l = 0), this becomes $-i \hat{S_{0}} \mathbf {k} / k^2$. Note that since $|\hat{\mathbf {R}}(\mathbf {k})| \propto k^{-9}$ and $|\hat{U}(\mathbf {k})| \propto k^{-3}$, we only consider the Brillouin zones with n ∈ [ − 1, 1]³ in Equation (5).

The resulting potential is computed using Equations (2) and (5), and then differenced to obtain the acceleration in Fourier space. An inverse FFT is used to obtain the mesh force in real space. Finally, the total acceleration field is computed by interpolating the acceleration field from the parent grid onto the refined grid, and adding it to the fine mesh force. The FFTs used for the refined meshes require zero-padding, but because the short-range force is compact, we have found that we only require an extra $0.65 \mathcal {R} a/\delta$ cells, where $\mathcal {R}$ is the usual refinement factor between levels. In this paper we always use $\mathcal {R} = 2$.

For illustrative purposes a decomposition of the gravitational force is shown in Figure 1 for a particle located at the center of the computational domain. The resolution of the root grid is 32³ zones, and we add three levels of refinement with respective coordinates [0.1875; 0.8125]³, [0.3125; 0.6875]³, and [0.40625; 0.59375]³. The components of the force on each level are shown in Figure 1, and compared to the analytical force F_anal∝r⁻², where r is the distance from the central particle. We can clearly see how each subgrid contributes on the smallest scales only, while the larger scale components are computed on the levels above. Note that the computed acceleration field is always cell-centered. It is therefore interpolated linearly when the Euler equations are solved in case a staggered mesh hydro scheme is used.

Zoom In Zoom Out Reset image size

Our implementation within the Enzo AMR code includes parallelization through the MPI library. This is relatively straightforward, as most of the density construction was already parallelized, with the largest change being communication of the coarse accelerations from the parent grid to the refine grid. In addition to parallelization, we have implemented a method to cache frequently used Green's functions, so that they do not need to be recomputed for each new grid.

We first verify that there is no self force acting on a particle. This is done by modeling the evolution of an isolated particle that is given an initial velocity and goes through the different levels of refinement. If the particle does not experience any self force, its velocity should not change.

The root grid resolution is 16³ zones and there are three static levels of refinement with coordinates [0.1875; 0.8125]³, [0.3125; 0.6875]³, and [0.40625; 0.59375]³. The particle of mass m = 1 is initially located on level 0 at coordinates r₀ = (0.15, 0.15, 0.15), and is given a velocity v₀ = (1.0, 1.0, 1.0). The system is evolved until t = 0.4, and we investigate two cases.

In the first case, we do not allow any subcycling and all levels are advanced at a constant timestep dt = 3 × 10⁻³. With the APM solver, the particle acceleration is found to be of the order 10⁻¹⁴ on the deepest level, therefore completely negligible.

In the second case, the system is evolved at a fixed timestep dt₀ = 2.5 10⁻² on level 0, and with subcycles on the refined levels. Using the added criterion for timestepping (Equation (10)), the value of dt₀ leads to the deepest level (level 3) being evolved with a timestep similar to the constant timestep chosen in the first case. We show in Figure 2 the error on the particle velocity Δv = |v(t) − v₀|/|v₀| for both solvers. The error for both methods is still small, of the order of 10⁻⁴, and the APM solver is slightly more accurate than the multigrid solver. However, allowing subcycles introduced an error due to the time-integration, error that is symmetric in all directions. We will witness this effect in Section 3.5 as well. We thus conclude that the APM solver does not directly introduce any self-force but that evolving the different levels with different timestepping decreases the time-integration accuracy. This is essentially because with subcycling, the force contribution from the different levels are computed at different times. Figure 2 also demonstrates that the multigrid solver does introduce self-forces, in particular for particles located near the boundary of a refined region.

Zoom In Zoom Out Reset image size

In this test, we compute the acceleration of nearly massless particles in a gravitational field created by a heavy central particle. The dimensions of the root grid are 32³, and we add one static refined grid. We set up at the center a particle of density ρ = 1.0, and n = 5 000 test particles of mass m = 10⁻¹⁰ distributed randomly in log r, where r is the distance from the central particle. We study two cases: one in which the refined grid covers the region [0.4375; 0.5625]³ such that the particle is at the center of the subgrid, and one in which the refined region is [0.5; 0.5625]³ such that the particle is deposited on the corner of the subgrid. In both cases, we evolve the system for Δt = 10⁻¹⁰.

The different force components are plotted in Figure 3. The tangential component of the force F_tan should be zero while the radial component F_rad should follow the analytical result F_anal∝r⁻², but is softened for radii about two cell lengths. The largest inaccuracy in the force calculation is reached at the boundary of the nested grid (r ≈ 2). Both solvers give an overall good result, as the relative mean errors of the radial and tangential accelerations are only of the order of a few percent. However, the transition is much smoother with the APM solver than with the multigrid solver. With the APM solver the overall noise in the radial force is reduced in comparison with the multigrid solver (Figure 3).

Zoom In Zoom Out Reset image size

Note that in the case where the massive particles is located at the edge of the refined grid, one can see how the smoothing of the potential starts at larger distances on the root grid than on the refined grid. Moreover, the acceleration of the central particle in the first case is ≈10⁻⁸ and ≈10⁻¹¹ for the multigrid and the APM solvers, respectively. In the second case, the respective accelerations are ≈10⁻³ and ≈10⁻¹⁰. We therefore confirm that there is no self-force introduced by the APM solver, as demonstrated in Section 3.1.

We investigate here the error in a spherical distributions for which we can analytically calculate the potential. These are the sphere of constant density (i = 1), the isothermal sphere (i = 2) and the Plummer sphere (i = 3). We note M and R are the total mass and radius of the sphere, respectively. The density distributions for the three cases are given by

$$\begin{equation} \rho _i(r) = \left\lbrace \begin{array}{@{} ll}\rho _0 & r \le R, \ i = 1 \\ \rho _0 \left(\frac{R}{r}\right)^2 & r \le R, \ i = 2\\ \frac{\rho _0}{\left[1 + (r/R)^2\right]^{5/2}} & r \le R, \ i = 3\\ \rho _{{\rm out}} & r > R, \end{array} \right. \end{equation} \tag{ 11 }$$

where ρ₀ = 1 and ρ_out = 10⁻¹⁰. The corresponding potentials are obtained using Gauss's law:

$$\begin{equation} \Phi _i(r) = \left\lbrace \begin{array}{@{} ll}-\frac{M}{2R} (3-r^2/R^2) & r \le R, \ i = 1 \\ -\frac{M}{R} [1-\log {(r/R)}] & r \le R, \ i = 2\\ -\frac{M}{R} \left(\frac{2\sqrt{2}}{\sqrt{1+r^2/R^2}}-1 \right) & r \le R, \ i = 3\\ -\frac{M}{r} & r \ge R. \end{array} \right. \end{equation} \tag{ 12 }$$

In each case, the dimensions of the root grid are 32³ and R = 0.3. We allow up to two levels of refinement and the domain is refined using the parameter MinimumMassForRefinement which we set to 10⁻⁶, 10⁻⁵, 3 × 10⁻⁶ for i = 1, 2, 3. Refinement occurs if the mass in a cell exceeds this parameter. The resulting hierarchies are plotted in Figure 4.

Zoom In Zoom Out Reset image size

The two solvers produce comparable results, therefore we show in Figure 5 the results for the APM solver only. In all three cases, the accuracy in both the radial and the tangential directions is at least at the percent level, except for the isothermal sphere at small radii. This is due to the fact that the gravitational acceleration is smoothed near the center because it diverges.

Zoom In Zoom Out Reset image size

In order to compare the gravity solvers with periodic boundary conditions, we study the case where the gas density is a sinusoidal function:

$$\begin{equation} \rho _{\rm SW}(x,y,z) = 2 + \sin {\left(\frac{2\pi x}{P}\right)}, \end{equation} \tag{ 13 }$$

which according to Equation (1) leads to the periodic potential:

$$\begin{equation} \Phi _{\rm SW}(x) = -4\pi \left(\frac{P}{2\pi }\right)^2 \sin \left(\frac{2\pi x}{P}\right), \end{equation} \tag{ 14 }$$

where P is the period. We use periodic boundary conditions and two static levels of refinement with coordinates [0.34375; 0.65625]³ (level 1), and [0.421875; 0.578125]³ (level 2). We perform simulations for P = 0.2 and P = 1.0, and with 64³ zones on the root grid. Solving the Poisson equation with periodic boundary conditions requires that there is no net-mass in the computational domain. We thus subtract the average density of the system (ρ_av = 2.0) from the GravitatingMassField before solving Equation (2).

The multigrid and the APM solvers give comparable results on this test. In Figure 6, we therefore show the computed force and the relative error on the force along on the x axis (Equation (14)) for the APM solver only. For the case P = 1.0, the APM solver is quite accurate because the density distribution varies on larger scales. For the case P = 0.2, the accuracy slightly decreases. This behavior is due to the combined facts that the density varies on smaller scales and that the acceleration is smoothed by the TSC deposition and the S(r) function, leading to an effective loss of resolution.

Zoom In Zoom Out Reset image size

Finally, we perform a two-body problem with the particles initially in a circular orbit in the x–y plane. The central particle has a mass M₁ = 1.0 and is located initially at the center of the domain. The test particle has a mass M₂ = 0.1 and is placed at a distance a = 0.3 from the central particle, at coordinates r₂ = (0.2, 0.5, 0.5). The resolution of the root grid is 16³ zones and we study four different cases.

3.5.1. Initial Configurations

3.5.2. Results

The trajectories of the particles are showed for the different cases in Figure 7.

Zoom In Zoom Out Reset image size

In all cases, the APM solver gives much more accurate results than the multigrid solver. With the multigrid solver, particles have already left the computational domain after less then five orbits (Figure 7, left column). Note that the particles follow periodic boundary conditions. One should also highlight that subcycling has very little effect on the multigrid solver.

On the contrary the APM solver yields a much more accurate evolution of the system (Figure 7, right column). Simulations in which subcycling is allowed (cases 1 and 3) show a small resonance of the system: the center of mass of the system is shifted during the evolution. This behavior is most prominent in the cases where both particles are located in the same refined level (Figure 7, right column, row 3). Note that while the particles are located in the same level of refinement, they are usually not covered by the same grid. This resonance is due to the fact that the level where the particles live (level 2) is evolved through 4 subcycles while the root grid level only goes through one, leading to a time-integration inaccuracy. More quantitatively, the orbital separation of the system has changed by about 3% at the end of the simulation (case 3). The center of mass has also moved by about 1% in the x and y directions.

If subcycling is turned off and all levels are evolved with the same timestep, the results become much more accurate with the APM solver (Figure 7, right column, rows 2 and 4). Regarding case 4, the orbital separation of the system has changed by 0.7% only. Moreover, the center of mass is stable within a relative error ≈10⁻⁵. We thus conclude that the force calculation with the APM solver is quite accurate, and that resonance in the orbit can be avoided by forcing the different levels to be evolved with the same timestepping.