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

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.

• As a remark to approach 1. You can mostly avoid cache misses by allocating the object data sequentially in memory (which is not easy but can be done). However you need to use general functions that handle all element orders. If you also want to include different types of elements you would have to deal with virtual dispatch but depending on the amount of work that you do per element this shouldn't be too expensive, and it should have a smaller effect as the element order is increased. – gnzlbg Aug 28 '12 at 13:24

## 3 Answers

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).

• That's what I thought as well. But will this efficiency still hold up if I need to get further details about the cells before I can process them? I mean, here I understand that I would at least $N_p$ for each cell to be processed, but there could also be a state variable (e.g. is_active_cell) that will affect how the cell is processed. – Michael Schlottke-Lakemper Aug 26 '12 at 6:52
• You would store each attribute you need to query in long arrays that you traverse sequentially. – Wolfgang Bangerth Aug 26 '12 at 12:53
• OK, but just to clarify: The long array with attributes is different from the array with the simulation data, correct? – Michael Schlottke-Lakemper Aug 26 '12 at 13:58
• @gnzlbg: What I meant is that if, say, your attribute is whether a cell is enabled or disabled (or which part of the domain it belongs to, or ...) requires us to store a bool (or a short int, or ...). On the other hand, the elements of the arrays we are talking here are, for example, the $N_p$ indices of the degrees of freedom on a given cell and are, thus, unsigned integers. Now, you can't store a bool in an array of unsigned integers -- it has to be a separate array. – Wolfgang Bangerth Aug 28 '12 at 12:43
• It is not necessary for data locality for everything to be in one array. It is quite ok from a processor/cache's viewpoint to have a small number of arrays that you traverse sequentially. The problems with having a "struct of everything that relates to a cell" and then an array of those are: (i) you group together variables of different data types that may need padding in between and aren't always accessed in the order in which they appear in the structure; (ii) in a given loop one often only accesses a subset of these attributes, leading to non-consecutive memory access and cache waste. – Wolfgang Bangerth Aug 29 '12 at 11:01

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).

• I wonder how exactly you handle random access to cells. For example, if you want to store the neighbor of a cell, how do you go about this? Do you use pointers? Or store two values ($N_p$ and id)? – Michael Schlottke-Lakemper Aug 28 '12 at 9:33
• @MichaelSchlottke I just store pointers to the neighbors. In my case (hierarchical quadtree/octree grids) it is not required because I can "cheaply" recompute the neighbors of a cell but I find its easier to store pointers in the cells. It's a lot of memory tho, about 4ptrs in 2D and 6ptrs in 3D. For exchange/window cells a ptr to the neighbor would not make sense (it lies in another memory space), but those cells have to be treated differently anyway. – gnzlbg Aug 28 '12 at 12:25
• I think that the scheme you describe in the answer will be approximately par for the course with option 2. I'd be surprised if they differed in any noticeable way. The only difference I see is that in your case writing loops over all cells is somewhat more awkward; since one writes a lot of loops over all cells, this means a lot of boiler-plate code. Other than that, this looks like a decent option. – Wolfgang Bangerth Aug 28 '12 at 12:41
• @WolfgangBangerth I haven't compared the performance of both approaches but I want to do that in the future. You are right in that the code required to do this does look like: "for_each cell<type1> in container<type1> do_stuff<type1>(..) ... for_each cell<typeN> in container<typeN> do_stuff<typeN>(..)". However, there is no boiler-plate because this is so repetitive that the compiler is able to generate it itself. DRY! – gnzlbg Aug 28 '12 at 13:19
• Could you show a small example in your answer of how this would work in C++, or is this asking too much ;)? – Michael Schlottke-Lakemper Aug 28 '12 at 13:46

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).