Data Structures and Algorithms

This page contains notes about the some of the complex data structures and algorithms in AREPO, which may be useful for development.

Voronoi Mesh

In Arepo (unlike in Gadget) the hydro domain is always between 0 and BOXSIZE in all dimensions. (Though the boxsize can be selected separately in each dimension.) The Voronoi mesh extends beyond this, though. The mesh is generated by inserting points into the initial tessellation, which consists of only one huge tetrahedron that is large enough that the hydro box is entirely contained inside it.

There are a few main data structures that the code uses, defined in voronoi.h.

DP - An array of points used for the mesh construction. (There are actually 5 points with negative indices, these are the vertices of the initial huge tetrahedron, and point -5 which is designated as an “infinity” point.) The point structure is not just a simple 3-vector, it has 2 3-double-vectors (x,y,z) and (xx,yy,zz), a 3-long long int-vector (ix,iy,iz), an id, and several more flags. The global variables Ndp denotes the current number of entries, and MaxNdp the allocated size of the array. As far as the function of the 3 fields, the meaning of the various coordinate entries in the DP structure is somewhat intricate. The coordinates are in fact repeated three times at the moment:

  1. x, y, z

  2. xx, yy, zz

  3. ix, iy, iz

The first triple, x/y/z, are just the ordinary coordinates, in the coordinate system and with the units that the simulation uses. The second triple, xx/yy/zz, are translated and scaled coordinates such that the floating point range of the resulting numbers falls into the interval [1.0,2.0[ for all points in the tessellation. Finally, ix,iy,iz are integer coordinates of the points, constructed by reading out the mantissa of the IEEE double precision numbers xx/yy/zz (this is done by the bit-fiddling operations). Because all xx/yy/zz have the same exponent (due to the 1.0-2.0 range) in the binary representation of the floating point numbers, these 52-bit mantissa integer numbers form a one-to-one mapping of all representable double precision numbers in the 1.0-2.0 range.

With the above construction in hand, AREPO has a robust and still reasonably fast way to deal with possible degeneracies in the Delaunay construction… the latter is actually a really major headache. The approach is basically the following: All the geometric predicates are evaluated first with the xx/yy/zz coordinates, using fast floating point arithmetic. Whenever there is a danger to get possibly the wrong sign in one of them, the code calculates the predicate with exact integer arithmetic based on the ix/iy/iz coordinates. Because the latter nicely map to the xx/yy/zz numbers, this is equivalent to evaluating the sign of the geometric predicate with arbitrary precions floating point arithmetic, except that the integer trick is faster and simpler to code. Note that because the mapping x/y/z -> xx/yy/zz can actually change the sign of some of the predicates, one cannot simply use the x/y/z here instead - perfect consistency can only be achieved between the xx/yy/zz and ix/iy/iz spaces. A downside of this approach is the storage required to hold multiple versions of the coordinates, but for the moment I have preferred this for speed reasons over a possible on-the-fly conversion of the coordinates whenever they are needed.

Note that the configuration option OPTIMIZE_MEMORY_USAGE reverses this speed <-> memory tradeoff, reducing the size of the point structure at the cost of additional arithmetic operations. It is generally not recommended except for specific circumstances.

DT - This is an array of Delaunay tetrahedrons (tetra structs). They contain the DP indices to the 4 points in field p, the DT indices of the neighboring tetras in field t (ordered such that t[0] is opposite to the point p[0]). The s field contains the p-indices (0-3) of the point in the neighboring tetra t[i] that is opposite to our point p[i]. That is, if we are sitting on tetrahedron i and want the opposite point in the tetrahedron that is opposite to our point j, we get it by DP[ DT[ DT[i].t[j] ].p[ DT[i].s[j] ]. The global variables Ndt denotes the current number of entries, and MaxNdt the allocated size of the array. The fact that the DT array stores the neighboring tetras means that you can walk your way through them. The functions get_tetra, InTetra, and image_get_next_tetra perform various variations of this.

DTC - This is an array of the centers of the “circumspheres” of the tetrahedra. These are the nodes where the faces of the Voronoi cells meet.

DTF - This is an array of char. Tetrahedron faces.

VF - An array of Voronoi faces (struct face). These contain two integers (p1, p2) which are the indices of the Voronoi cells on each side of the face, the face area, and (cx, cy, cz) the center coordinates of the face. The global variables Nvf denotes the current number of entries, and MaxNvf the allocated size of the array. The main hydro loop (flux computation) is done by looping over this array.

DC - This is an array of edge* which encodes the connectivity of the Voronoi cells as a linked list. Each entry corresponds to one edge in the directed graph of Voronoi cells and edges between them. It has indices index and dp_index, which indicate the into the SphP and DP arrays of the mesh point on the other side of the face. For a given cell in SphP, the first_connection member points to the first entry in the DC array corresponding to edges going out from that cell. The next member is the DC index of the next edge, and SphP.last_connection is the index of the last entry in the DC array corresponding to the cell. (Note that you can not look for the end of the linked list, because even the last_connection entry has a valid next member pointing to an edge corresponding to another cell. Why is this?)

image_flags

The code constructs mesh cells that overlap processor boundaries and the box boundary through ghost points. In particular, the boundary conditions at box borders are implemented by inserting primary mesh generating points (the ones inside the box) if needed several times by translating them periodically in the x/y/z directions (or mirroring them for reflective boudaries). Only one such translation in each dimension is allowed, so that there is effectively a 3x3x3 replication grid of the principal box, and each mesh generating point could be at most inserted 27 times into the mesh.

The bitflags are basically used to encode whether a principal mesh-generating point should be inserted as is into the mesh (unshifted), or whether it should be shifted by -Boxsize or +Boxsize in one or several spatial dimensions. Because the code needs to keep whether a given point has already been inserted into the mesh before, it converts the position value from above (possibilities 0-26) into a bit position, so that through one 27-bit wide bitfield stored in a 32-bit integer all possible combinations of multiple insertions in the 3x3x3 box mesh can be encoded.

As far as the DC[] objects are concerned, each principal point stores its neighbors separately, and whether they are unshifted or shifted. In the domain decomposition, it can now happen that a principal point is shifted by -BoxSize or +BoxSize in one (or several) dimension because the periodic wrap around wants to map it back to the principal domain [0,L]. It’s neighbours listed in DC should then also be shifted accordingly, which boils down to modifying their bitflags (assuming that they are not also wrapped, which may happen, but would not lead to problematic consequences) accordingly.

A Terminate(“a”); or Terminate(“b”); means the code wanted to alert us of an unexpected situation, where apparently a point that wants to be wrapped by +L (because it has moved to a negative coordinate value) happened to have a neighbor that was shifted by +L already from its principal box coordinate (and hence one would think that it should now be shifted by +2L after the wrap). However, it should be perfectly safe to simply ignore this situation by disabling the Terminate(); or replacing it with a warning. The next mesh construction will rectify the situation by finding new correct DC’s for the affected cell.

Domain decomposition and ghost cells

Due to the domain decomposition needed when running on more than one processor, the points must be partitioned, but points on the boundary need to be duplicated so each processor can generate a correct tessellation of the cells in its domain. This is indicated by the DP[].index and DP[].task fields of the point structure. There are three different types of Delaunay points:

  1. “local ones”

  2. “local ghosts”

  3. “foreign ghosts”

The DP in category (a) have DP[].task == ThisTask, and their coordinates are identical to the coordinates of a local mesh-generating point identified by DP[].index (i.e. P[DP[].index].Pos), with the corresponding hydro quantities given by SphP[DP[].index]. Such points must have 0 <= DP[].index < N_gas, and their corresponding Voronoi cells are “primary cells”. The sum of the volumes of all primary cells over all processors is the simulated volume.

Now, it can also happen that due to periodic boundaries or reflective boundaries, a local mesh-generating point needs to be inserted into the local tessellation a second (or third…) time, forming a local ghost point/cell. This will for example always happen if a point shares a wall with a reflective boundary. Such points are in category (b) and have DP[].task == ThisTask but DP[].index > N_gas. The hydro-quantities of the ghost cell can be accessed without communication in this case, since the primary cell is on the same MPI-task. My convention is that these hydro-quantities can be accessed through SphP[DP[].index - N_gas].

Finally, some Delaunay points are ghost points whose origin lies on different processors. Here DP[].task simply gives the task-rank of this foreign CPU, and DP[].index is now used as an index into the local PrimExch array, where the data of primary Voronoi cells from the foreign CPUs is imported.

There are two more fields in the DP array related to domain decomposition. The originalindex and ID entries in DP[] were introduced to change how fluxes across Voronoi faces are treated that appear in the local tessellations of two (or several processors). If a cell i is adjacent to a cell j, but the corresponding mesh-generating points lie on different processors A and B, then both processors will construct the face between i and j as part of their local Voronoi tessellation. They will also mutually import the other point as a foreigh ghost point. The question then becomes how the flux across the face should be calculated such that a possible double counting or the like is avoided. Now, the code unambigously decides that one of the processors involved in the shared face should take responsibility for it. This is accomplished through the ID of the Delaunay points (and this is also one reason why only unique IDs should be used for all points). Only this responsible CPU calculates the flux for the face. But it is then necessary to “apply” the flux to a foreign ghost point if needed. The “originalindex” is the index of the SphP[] structure corresponding to an imported ghost point on a foreign processor. Knowing this, one can tell a foreign CPU which of its cells should receive a certain update.