(How) do you take into account memory fragmentation?

I use an example from finite element theory, but anybody who maintains a large datastructure and successively extends it will find something similar.

Suppose I have an unstructured mesh of points and triangles, where the points are given by coordinates (say $x$ and $y$) and the triangles each consist of three point indices (say $i$,$j$ and $k$).

As common in FEM, the mesh will be successively refined. If we resort to global regular refinement, the number of triangles will grow by factor $4$ with each iteration of the refinement. Depending on how this is done, the memory layout will develop differently.

Say the mesh occopies memory cells 1 to 300, anything beyond there being free.

Example 1:

We allocate the space for the new mesh, cells 301 to 1501, fill it up with the data of the refined mesh, and forget the old one. The next refined mesh will be placed in cells 1501 to 6300, the next one in 6301 to 21500, and so on. The location of the current mesh will move in memory "to the right", while a huge patch will not be used. We may run out of memory prematurely.

One might observe in the above example, that this will only hinder us for one refinestep, because even without that fragmentation we would run out of total memory one refinement later. As the vertex array is taken into account, too, the problem can become more severe.

How can this be circumvented?

Example 2:

Realloc the triangle array to cells 1..1200. Create the new mesh in cells 1201 to 2400. Copy the content of that working copy to cells 1..1200, and forget the working copy. Repeat similarly.

Ok, still we run out of memory prematurely, because we need a working copy. How about this:

Example 3:

Realloc the triangle array to cells 1..1500. Copy old mesh to 1201 .. 1500. Create new mesh in cells 1..1200. Then forget the copy of the old mesh.

The case here is artificial, because one would not use global mesh refinement on these scales. If the growth is much smaller, memory realignment is possible to avoid fragmentation. However,

Questions:

1. Does memory fragmentation ever become critical in practical scientific computing/high-performance computing?

2. If at all, how do you avoid it? Maybe my machine model is even wrong, and the OS by some heavy magic does realign the memory tacitly, or manages fragmented blocks on the heap.

3. More specific, how does it impact grid management?

I disagree with Matt about total memory use for the mesh and about extra indirection. For explicit methods, it is common for a very small number of vectors (e.g. 2) to represent all simulation state. For a scalar problem, just defining the coordinates of the mesh may be more than this, and if connectivity is explicit, it is substantially more. High-resolution simulations with explicit methods frequently need a huge number of time steps, so it's common to use reasonably small subdomains just to be able to complete a model run in a reasonable amount of wall time. Thus, it's not all that common to run out of memory because of the cost of storing the mesh.

A far more fundamental problem is the memory bandwidth requirement for each step of an algorithm (e.g. evaluating a residual or taking a time step). For unstructured methods, this involves some mesh-based information, but it usually does not require accessing all mesh storage. Note that even if solver memory dominates mesh memory (as is common for many implicit methods), the simulation is still frequently dependent on the bandwidth requirements of a mesh traversal. There are two factors:

What data is needed during routine mesh traversal?

In finite element methods, an element-to-global map is needed. For higher than first order, this is often implemented by associating degrees of freedom with intermediate entities like faces and edges, but once the element-to-global map is constructed, the intermediate entities are not needed. In particular, the connectivity of the intermediate entities is never used directly by a simulation. That information might not be touched for many iterations of an implicit solver or time integrator, but still be required during adaptivity or to set up a new function space. It is common to set up the element-to-global map in a simple array and not touch the "mesh" data structure during this period, in which case the storage format and data locality of the mesh itself is less important.

Finite volume methods require the cell volumes, face-to-cell connectivity, face area, and face normals. Note that vertex coordinates or connectivity need not be available. In extreme examples (e.g. FUN3D), all other information is discarded during offline preprocessing (mesh generation) and only this reduced representation is available to the simulation. This is very efficient, but precludes moving meshes and adaptive refinement.

How is that data laid out in memory?

For many simulations, with modest orders of accuracy and complexity of physics, performance is limited by memory bandwidth. Modern CPU architectures from IBM, Intel, AMD, and NVidia support an arithmetic intensity between 4 and 8 flops/byte. With such efficient floating point units, we should try to optimize our algorithms for memory bandwidth. This might involve somewhat deep algorithmic changes such as favoring unassembled high order methods (compare the "assembled" and (unassembled) "tensor" lines in the second figure), but we can start by attempting to fully utilize the bandwidth provided by the hardware. This often involves computing orderings (of vertices, faces, cells, etc) such that cache is reused as well as possible. It also involves exploiting blocking and ordering unknowns so as to have minimal metadata and to activate the right number of data streams. (Usually 3D unstructured simulation ends up with "too many" streams, but some modern hardware like POWER7 needs many streams to saturate bandwidth, so it occasionally makes sense to intentionally organize data to activate more streams.) The PETSc-FUN3D work provides a classical, but still highly relevant discussion of these performance optimization for unstructured implicit CFD.

1. Do not be afraid of malloc(). There is no need to pack the current refinement of the mesh into the same array as the coarser inactive part of the mesh. Whenever you no longer need an old part of the mesh, just free() it.

2. After refinement, compute a good ordering for that refinement level (Reverse Cuthill-McKee is popular, but more cache-specific orderings can also be used). If your load imbalance is reasonable (less than 10%, say), you can compute this new ordering locally (without parallel partitioning and redistribution). The cost will likely be similar to a single "physics" mesh traversal, but can speed up those traversals by a factor of 2 or more. This step will usually pay off any time mesh adaptivity does not occur on each "physics" mesh traversal. If such frequent adaptivity is required for your problem, you are likely making small changes and should still reorder occasionally. I would still avoid fine-grained memory pools because it makes reordering more difficult, but you can use reasonably large chunks to balance peak memory usage, ordering cost, and cost of incremental updates.

3. Depending on the discretization, consider extracting a reduced representation with low memory bandwidth requirements for the "physics" traversals. Unnecessary indirection is bad because it increases the bandwidth requirements and, if the target of the indirection is irregular, causes poor cache reuse and inhibits prefetch.

• This problem over intermediate connectivity taking up bandwidth is trivial to eliminate for your case, and is done is PETSc DMComplex. You divorce storage of mesh topology from field data. Storing element-to-global is wrong in every case I think. The topology has the same information in smaller space, and is much richer. You compose a traversal, like closure, with offsets into storage. Mar 3, 2012 at 16:38

In deal.II, we refine the mesh by throwing away old cells and replacing them with new ones. But new ones are also placed into the memory holes left by deleted cells. All loops over all cells are then done in the order in which cells are encountered in memory to keep cache hits high.

The bigger question is how you store the data that defines cells. You could of course do struct Simplex { Vertex vertices[4]; int material_id; int subdomain_id; bool used; void * user_data; };

class Triangulation { Simplex *cells; }; but this is not cache efficient because most loops over all cells will only touch a subset of the data you store in your Simplex data structure and so only a fraction of the data that lands in the cache will actually be used. A better strategy is to do something like this: class Triangulation { Vertex *vertices; int *material_ids; int *subdomain_ids; bool *used_flags; void* *user_data; }; Because in loops over all cells subsequent iterations are likely to access the same subset of data that defines cells, read-ahead caches will only preload that data that you're actually going to use and will consequently lead to high cache hit rates.

1) No. Solver memory far far far outweighs mesh memory. Even if you run the lightest, explicit solver, the mesh is at most 25% of the memory of the simulation, and much more likely < 10%.

2) Break down your allocation and use memory pooling. You do not need to allocate a contiguous chunk for the entire mesh since you usually only need to iterate over local pieces. Introducing one level of indirection does not impact the performance in a meaningful way.

• Any references for this claim? A simple matrix free method should have a very light solver footprint. Dec 20, 2011 at 18:46
• Sure, if you do matrix free, that may be a different issue. But short of that it's easy to see that Matt's comment is right: for each degree of freedom on the mesh, the matrix has to store as many doubles as this DoF couples with -- which in all non-trivial 3d cases runs into the hundreds (count this once for a Q2-Q1/Taylor-Hood Stokes element in 3d). You'll easily see that this takes much more memory than the data that defines the mesh. Dec 21, 2011 at 19:52

If you're running out of memory, just run on more nodes so you have more memory. You're wasting a very valuable resource (the human brain) to solve a problem that has a very easy solution.