How to efficiently structure simulation data in memory for cells with varying degrees of freedom?

Question

For a discontinuous Galerkin-based simulation I need to store cell-based simulation data in memory. Since the order of the polynomial approximation $N_p$ may vary between cells, I wonder what the most efficient data structure would be, if the typical access pattern is a loop over all cells? These are the ideas I had so far:

1) Use an array of cell objects, each object with a pointer to its (individually allocated) simulation data array
Pro: Very easy to handle algorithmically, no headaches when pre-allocating memory for data storage, very easy to store topological relationships between cells (i.e. neighbor information), straightforward direct access of specific cells.
Con: Very (extremely?) slow when looping over all cells as there will probably cache misses every single time.

2) Use an array of cell objects, each object with a pointer to its simulation data array (memory pre-allocated and cell data stored contiguously)
Pro: Same as above. Also probably faster since data is now contiguous in memory.
Con: Still slower than just using arrays, since each cell access still requires the knowledge of $N_p$.

3) Use an array (or an STL container) for each set of cells with common $N_p$
Pro: Each container can contain the data contiguously, thus there will be many cache hits if iterating over $N_p$ first, cell second.
Con: Difficult if memory pre-allocation is desired. Random cell access using only the cell id becomes very expensive.

Does anyone have experience with one of the aforementioned methods, or has an even better alternative in mind? If you need more information on e.g. the data access patterns, please let me know.

Note: This question is similar to this one on Stackoverflow, but rephrased and simplified.

Answer 1

deal.II uses variants of your option (2). This is of course also essentially the way you store sparse matrices and is very efficient since you only have two arrays (the array that for each cell stores the beginning index of the data that corresponds to this cell, and the array that contiguously stores all the data of all cells). You traverse both arrays sequentially, producing few cache misses. Both of your options are either terribly inefficient (option 1) or rather wasteful with non-consecutive access (option 3).

Answer 2

I use a variant of option 3, that is: N independent containers, one for each cell/element type/order such that all cells of the same type are aligned in memory. When I want to perform operation x on all cells (integrate, quadrature, evolve, numerical flux...) I do:

for each container<cell_type> in cell_containers
  for each cell in container<cell_type>
     do x<cell_type>(cell)

where x is a different, specialized and optimized function for each cell type. To allow random access a small vector of pointers to all cell containers is stored, inside each container random access iterators are available (which should be at least as fast as the second approach but here your jumps in memory are constant inside each container because the elements of the same type are aligned sequentially). If you need global ordering for say a space-filling curve I can't see how implementing it to iterate over elements[order][id] should be more difficult than elements[id] but I haven't tried it.

The main advantage is that the compiler knows a lot about what is going on and can safely optimize (no virtual dispatch, no aliasing, no breaks out of the loop, the types and element order inside your numerical kernels are all compile-time constants...). As a nice side effect, hybridizing an MPI code using this approach is just as easy as loading the std parallel library(*) and compiling with OpenMP, no extra code required. For boundary elements one has usually to perform other operations and store different data. I store them separately using this same approach and define specialized kernels for them.

The main drawback is that you need to perform memory allocation very carefully. The problem is, before hand the only thing you know is how much memory is available, but you don't know how many elements of which order you are going to need. In particular, imagine that you have allocated a certain amount of memory for a certain type of cell, and you want to exceed it. This will require that you have at least enough free memory to copy the whole vector of cells plus one extra cell. It also involves copying all the cells from that vector to the new one but that is fast. This is not smart, you probably want to allocate memory for more than one extra cell, but for how many? You probably dont want to allocate memory for twice as many cells. If you want to do this in a smart perfomant way you really need your own allocator that knows about your data structures, your system, and the on going simulation.

Pros:

• Operations that need to be performed on all cells (numerical kernels) are very fast.
• Random access is available and is as fast as with the "one vector of cell data/one vector of cell ptrs" approach.

Cons:

• Memory allocation has to be performed very carefully (but isn't this always the case?).

_{* using TBB requires replacing your for_each loops with parallel_for_each or you can just implement your own for_each that takes a parallelization policy (or macro it: please dont do this).}

Answer 3

Also think about the loops over faces where you compute the flux once and update the left and right cell (contributions LL,LR, RL and RR in the matrix).