Neighbor Lists

The sums used in the SPH expressions above are dependent on each particle's kernel function, which in turn depends on each particle's smoothing length. How do we determine what it should be, and how it will evolve over time? Our method works by defining an optimal number of neighbors for each particle. This value is set by the user in the input files (see Sec. 5.1.1) by the parameter NNOPT. The kernel function of the particle, a sphere with radius equal to 2 smoothing lengths, ideally will overlap precisely this number of other particles. In practice, we use a relaxation technique, such that a every iteration of the code, the smoothing length is set equal to the arithmetic mean of the previous value and the power law estimate of the new value, thus

$$h_{new}=\frac{1}{2}h_{old}\times \left[1+\left(\frac{N_{opt}}{N_N}\right)^{\frac{1}{3}}\right].$$ (9)

In many cases, it is advisable to place limits on the minimum and maximum smoothing length allowed in the calculation. In practice, we have found that a minimum value is generally unnecessary, and setting a fixed floor could lead to particles in high-density regions reaching the maximum number of neighbors allowed in memory before all overlapping particles are reached. Setting a maximum is often important, since smoothing lengths grow dramatically when particles find themselves alone, and in this limit the SPH density kernel inherently does a poor job of describing the local density field. Typically, setting the maximum smoothing length to be comparable to the stellar radius is appropriate. The maximum and minimum smoothing lengths are set with the parameters ``HMAX'' and ``HMIN'', respectively. Setting ${\tt HMIN}<0$ tells the code there is no fixed minimum for the smoothing length.

To determine which particles neighbor each other, we use a linked list/head-of-cell method. First, we define both an integer vector with one entry for every particle and a 3-dimensional integer grid around the matter, whose size is set by the parameters NCMX, NCMY, and NCMZ in the file nelist.f. The physical extent represented by each grid cell in each direction is calculated by taking the average smoothing length for all particles, and multiplying by a factor of $1.3$ (this just works faster). Once the grid dimensions are defined, we run a loop over all the particles in the simulation, calculating the grid cell in which each one lies. If the particle is the first one to fall in its grid cell, we enter its number into the grid. If there are already one or more particles in the grid cell, we copy the old particle's number to the vector in the new particle's slot, and enter the new particle's number in the grid cell. When we are done, the algorithm to find every particle in a given grid cell is simple: the number written in the grid cell, $i_1$, is the highest-numbered particle which lies in the cell. The $i_1$th element of the vector contains $i_2$, the next-highest numbered particle in the cell. By making our way through the vector, we recover all particles in the cell, until we reach an element in the vector containing zero, indicating that the cell has been accounted for completely.

Next, we go through the process of assembling the neighbor lists for each particle. For each particle, we loop over all grid cells which may lie within two smoothing lengths of the particle, and determine for every particle in the cell whether this is the case. For each particle meeting the criterion, we record its value in the array ``NNI''. Since all hydrodynamical interactions are parallelized as described in Sec. 3.1.3, we only record the neighbor lists spanning particles with ${\tt n\_lower} < i < {\tt n\_upper}$ on each processor. We also keep track of the number of neighbors each particle has, ${\tt NN}$, but in a way such that every processor has access to the complete list for all particles.