Unstructured sparse dimension

Neighbor chain

In order to represent locations which are indirect neighbors of a starting location, that is neighbors of neighbors (of neighbors …) of the starting cell/edge/vertex, we use a list such as edge->cell->vertex or cell->edge->cell->edge->cell. The direct neighbors are simply edge->cell, cell->vertex and so on. The same element twice in a row, e.g. cell->cell is not allowed, since it is simply a shorthand (for cell->edge->cell in this case) and hence redundant. The target location type is the last element of the list. All locations of the target type traversed while following the chain’s links are selected as indirect neighbors, excluding the starting location (which may or may not be of the target type). See figure below.

A partial ordering when returning the neighbors is adapted by any mesh interface compliant with dawn. The key of the ordering is simply the distance d(X,Y) in the graph defined by the mesh. See the following examples:

The Diamond:

E->C->V = {{v₁, v₂}, {v₃,v₄}} since d(e,v₁) = 1, d(e,v₂) = 1 and d(e,v₃) = 2, d(e,v₄) = 2 (since the green vertices v₃, v₄ need an additional "hop" via the cells, whereas the purple vertices v₁, v₂ are connected to edge e "directly".

V->C->E = {{e₁, e₂, e₃, e₄, e₅, e₆}. {e₇, e₈, e₉, e₁₀, e₁₁, e₁₂}

Sparse dimension

Needed to store data, which is local to a starting location, on the (direct or indirect) neighbors of that location.

A field can be defined on an unstructured dimension (UnstructuredFieldDimension). A dimension would have a dense part and an optional sparse part. By abuse of notation we call a dimension with sparse part specified a sparse dimension. The dense part tells on which location type the field is defined. If no sparse part is specified, the field will have one value for each dense location. The sparse part is a neighbor chain whose starting element is the same location type as the dense part’s. When specified, the field will have, for each dense location, as many values as locations selected by the sparse part. This corresponds to storing, locally to a dense location, an array on its direct/indirect neighbors (specified through the neighbor chain).

In a neighbor reduce, we would accept as rhs either a field with only dense part ( = reduce over direct neighbors), or same dense part as lhs’s type but non-empty sparse part ( = reduce over indirect neighbors defined by sparse part).

To write into a field with non-empty sparse part, we can use the assignment expression, but with constraints on the rhs: the rhs must be an expression whose fields accessed (if any) have either dense part same as lhs’s target type (last location type of sparse part) and no sparse part or dense part same as lhs’s starting type (first location type of sparse part) and no sparse part or exactly same dense and sparse parts as lhs’s. Such assignment expression translates into a loop (over indirect neighbors) nested into the parallelizable loop (“over the dense”). (See code generation examples at the end)

Using sparse dimensions in combination with weights makes it necessary (at code generation) to know the dimensionality of neighbor chains beforehand (in order to generate meaningful error messages). This is only possible in structured meshes. For a hexagonal triangular mesh the mapping would be:

|cell->vertex| = 3
|cell->edge| = 3
|vertex->edge| = 6
|vertex->cell| = 6
|edge->vertex| = 2
|edge->cell| = 2

Currently, no such checks are done at the point of writing (neither at compile nor runtime) and the user is responsible to put the correct number of weights.

Examples

E1:

gridtype triangular;
field(e->c) a;
field(e) b, y;
field(c) z;
...
a = b * z;
y = reduce(e->c, op=+, init=0, expr=(a), weights=[1, -1]);
...

----- gets code generated into -----

const W[2] = [1, -1];
for(e : edges) {
    int i = 0;
    for(c : cell_neighbors(e)) {
         a(e, i) = b(e) * z(c);
        ++i;
}
int i = 0;
y(e) = 0;
for(c : cell_neighbors(e)) {
         y(e) += a(e, i) * W[i];
        ++i;
}
}

E2:

gridtype triangular;
field(e->c->v) a, b;
field(e) x, y;
field(v) p;
...
a = b - p + x + 2;
y = reduce(e->c->v, op=+, init=0, expr=(a * (b + p)), weights=[0.5, -1, -0.5, 1]);
...

----- gets code generated into -----

const W[4] = [0.5, -1, -0.5, 1];
for(e : edges) {
    int i = 0;
    for(v : indirect_neighbors(e, e->c->v)) { // neighbor chain
        a(e, i) = b(e, i) - p(v) + x(e) + 2;
        ++i;
}
int i = 0;
y(e) = 0;
for(v : indirect_neighbors(e, e->c->v)) { // neighbor chain
    y(e) += (a(e, i) * (b(e, i) + p(v))) * W[i];
    ++i;
}
}