Set Algebra and Recursive Interval Representation

A Unified Framework for MultiDimensional Adaptive Mesh Refinement Methods

Loïc Gouarin

HPC@Maths team
CMAP / CNRS - Ecole polytechnique

Context

Patched-based representation

Advantages

• compact data structure
• rectangular zones work well for tiling and optimizing caches
• efficient to work with neighbors and the different levels

Disadvantages

• requires more cells than necessary
  • grid hierarchy: overlapping
  • larger refinement zone

Cell-based representation

Advantages

• no grids hierarchy
• keep only needed cells

Disadvantages

• more difficult to find the neighbors
• more difficult to work with the different levels
• tree-like structure
• more difficult to add ghost cells hierarchy

If we look at all the open source software specializing in dynamic mesh adaptation, there are two main families: - patch-based, which is a hierarchical representation of the mesh: layers are placed on top of layers - cell-based, which is a flat representation of the mesh

Each has its advantages and disadvantages: - patch-based has rectangular zones for tiling and optimizing caches. But it generally requires more cells than necessary. - cell-based requires far fewer cells but requires a tree-like structure, which means that you lose the good memory locality you had with patch-based. We use space filling curves such as Morton or Hilbert to find an acceptable locality.

Adaptive methods

AMR

Advantages

• based on heuristic criteria (gradient, second derivative, …)
• easy to implement

Disadvantages

• physical problem understanding is important
• add potentially more cells that needed

MRA (see the presentation of M. Massot)

Advantages

• based on wavelet decomposition
• error control
• independent of the physical problem
• only needed cells

Disadvantages

• difficult to implement with the current data structures
• can be costly

Can we have a data structure that takes the advantages of both and well suited to implement various adaptive mesh refinement methods ?

Samurai

Design principles

• Compress the mesh according to the level-wise spatial connectivity along each Cartesian axis.
• Achieve fast look-up for a cell into the structure, especially for parents and neighbors.
• Maximize the memory contiguity of the stored data to enable caching and vectorization.
• Manage and facilitate complex operations between meshes (numerical schemes, MRA, …).

An overview of the data structure

[ start , end [ @ offset

We illustrate the samurai data structure, starting with a 2D problem to see all the components. The data structure is fully recursive and can therefore handle Nd Cartesian meshes.

The main element of samurai is an interval modeled by its beginning and non-included end. There’s also an offset parameter that will serve us in two different ways. We’ll come back to this later.

An overview of the data structure

• Level 0
• x: [0, 4[, [0, 1[, [3, 4[, [0, 1[, [3, 4[, [0, 3[
• y: [0, 4[
• y-offset: [0, 1, 3, 5, 6]
• Level 1
• x: [2, 6[, [2, 6[, [2, 4[, [5, 6[, [2, 6[, [6, 8[, [6, 7[
• y: [2, 8[
• y-offset = [0, 1, 2, 4, 5, 6, 7]
• Level 2
• x: [8, 10[, [8, 10[, [14, 16[, [14, 16[
• y: [8, 10[, [14, 16[
• y-offset: [0, 1, 2, 3, 4]

The first thing we’re going to do is go through each level and each y-component, and look at the set of contiguous cells along the x-axis to construct intervals.

Identify the different types of cells

Identify the different types of cells

Identify the different types of cells

Identify the different types of cells

Algebra of sets

\

=

Algebra of sets

The search of an admissible set is recursive. The algorithm starts from the last dimension (y in 2d, z in 3d,…).

The available operators in samurai are for now

• the intersection of sets,
• the union of sets,
• the difference between two sets,
• the translation of a set,
• the projection of a set.

The search for an admissible set is completely recursive. We’ll start by searching in the largest dimension and if we don’t find anything we’ll stop.

If an admissible interval is found, we look at it in the lower dimension, which means we manipulate very few intervals because we’re only looking at a subset.

Usage example: find jumps

auto jump_set = intersection(
                    translate(mesh[1], {1}).on(0),
                    mesh[0]
                );

Usage example: the projection operator

auto proj_set = intersection(mesh[level], mesh[level + 1])
          .on(level);

proj_set([&](auto i)
{
    u(level, i) = 0.5*(u(level+1, 2*i) + u(level+1, 2*i+1));
});

Compression rates

Compression rates

Level	Num. of cells	p4est	samurai (leaves)	samurai (all)	ratio
\(9\)	66379	2.57 Mb	33.68 Kb	121 Kb	21.24
\(10\)	263767	10.25 Mb	66.64 Kb	236.8 Kb	43.28
\(11\)	1051747	40.96 Mb	132.36 Kb	467.24 Kb	87.66
\(12\)	4200559	163.75 Mb	263.6 Kb	927 Kb	176.64
\(13\)	16789627	654.86 Mb	525.9 Kb	1.85 Mb	353.98
\(14\)	67133575	2.61 Gb	1.05 Mb	3.68 Mb	709.24

We can see here that the mesh is very compressed compared to a version using a tree structure.

The total size of the mesh is multiplied by just under 4, as we build different meshes to handle ghost cells and prediction cells.

What’s important to note is that in a parallel context, it doesn’t cost much for each process to have the meshes of its neighboring domains.

Other features

• Loop algorithms over the levels and the cells
• Simplified access operator
• Reconstruct a value anywhere even if the cell doesn’t exist
• Helper classes to construct complex meshes
• Helper classes to construct schemes for explicit and implicit usage
• Helper classes to construct N-D operators and expressions using xtensor or Eigen
• MPI support
• HDF5 support

Roadmap

Here’s the samurai roadmap for the coming months and probably years.

We now have the possibility of changing the containers in which our fields are stored. We did this so that we could start using Kokkos. So we’d be very interested in working with the people at Cexa on this.

We have received several funding packages that have enabled us to hire two research engineers.

As part of numpex, we also have an engineer position to work on the I/O part of AMR methods in collaboration with Dyablo developers.

People involved

• 4 core developers
• users and developers communities
• Various collaborations and research programs

Real-World Applications

• Interfacial flow simulation
• Simulation analysis on the Hydrogen risk
• Plasma dynamics
• Battery simulation
• Combustion modeling

Samurai

https://github.com/hpc-maths/samurai