Mesh traits and concept helpers

[havogt]

• is_connectivity

• maybe has_connectivities (has_iterationspace)

• type for max_neighbors etc

• move to mesh.hpp and define LocationType tags (namespace topology)

[anstaf]

The design challenges:

• As for today the mesh concept assumes that any connectivity that is defined by neighbour chain can be queried from the concept model. This is convenient for the concept users but pretty painful for concept modelers.
• Another related problem: how to write is_mesh trait? What connectivities should we check? all possible pairs of locations? And additionally all possible triples (like cell-cell-edge, cell-cell-cell etc. )? How about quadruples? Should me explicitly ban them? Or support?
• Should we check that the mesh keeps fundamental topological invariants? For example that maximum neighbour number of the edge->vertex connectivity is 2. Or maximum neighbour for vertex->vertex equals the one for vertex->cell?
• How about reduced topologies? For example the mesh that supports only vertices and edges (but not cells) in its neighbour chains. The extreme case here: we have only vertices and they form cartesian grid. Is it useful use case?
• If the model of the mesh doesn't provide connectivity for the given neighbour chain, should the mesh::connectivity generate the default? It looks like that for the triples, quadruples etc. it is possible to generate neighbour table on the fly if all pairs n the chain are available. Even more: if there is for example cell-> edge and edge->vertex, cell->vertex is generated in linear time (against cell size).

[anstaf]

Observations:

• cell->cell == cell->edge->cell ( cell->vertex->cell also builds cell to cell connectivity. It is covensional what is the right one).
• vertex->vertex == vertex->edge->vertex

• edge -> edge == edge -> vertex -> edge

• cell -> vertex == cell->edge->vertex
• vertex -> cell == vertex -> edge -> cell
• max_neighbour(A -> B -> A) == max_neighbour(A -> B) * (max_neighbour(B -> A) - 1)
• max_neighbour(As ...-> B -> C) <= max_neighbour(As ...-> B) * max_neighbour(B -> C)
• max_neighbour(edge->vertex) == max_neighbour(edge->cell) == 2
• max_neighbour(vertex->edge) > 1
• max_neighbour(cell->edge) > 2
• max_neighbour(vertex->cell) == max_neighbour(vertex->edge)
• max_neighbour(cell->vertex) == max_neighbour(cell->edge)
• edge->vertex neighbor table can not contain skip values

[anstaf]

Let us restrict the requirements to the neighbor chain with the following:

• the left and the right hand side of every link in the chain must be of different locations.

I.e the chain vertex->vertex->vertex is illegal. vertex->edge->vertex->edge->vertex should be used instead (which is semantically the same)

[anstaf]

Should Euler formula for the mesh be validated somewhere in the concept helpers? It is clear that meshes on torus or Klein bottle are not the target of the concept. However we can distinguish between meshes on sphere and meshes on disk (where the "outer" cell is explicitly excluded from neighbour tables). Should we support both?

[anstaf]

Let us assume that we calculate something on the mesh on disk. Suppose we are taking connectivity for the chain that ends on cell. For example vertex->edge->cell. If for the given vertex number the neighbours list contains -1 (skip value), the semantic meaning of that is ambiquis. Either the topological number of the neighbours for that vertix is less than max_neighbour, or there is a neighbour that on the other side of the disk border.

[anstaf]

Mesh Contracts

Mesh API for generated code

Generated code accepts mesh as a template parameter.

template <class Mesh, class... Fields>
void nabla(Mesh, Fields&&... fields)  {...}

Here Mesh is an empty class. This means that it is cheap to pass it everywhere by value (including kernel functions). The "real" state of the mesh is hold as a singleton. This singleton stays immutable after initialization.

Things that can be performed with the instance of Mesh (aka obj):

• mesh::init(obj) - initialization of the internal singleton state.
• mesh::get_location_size<Location>(obj) - returns the number of topological locations as integral constant
• mesh::get_neighbor_table(obj) - returns a SID with the neighbor connectivity table. The supported syntax:

  mesh::get_neighbor_table(obj).from<cell>().to<edge>();
  mesh::get_neighbor_table(obj).from<cell>().to<vertex>().via<edge, vertex, edge>();

• mesh::field_builder(obj) - create a field. The syntax is gridtools::storage-like:

  mesh::field_builder(obj)
    .type<double>()
    .at<vertex>()
    .extra_dimensions<dim::k, dim::another>(80, 4_c)
    .initializer([](auto v, auto k, auto another) {return v * 42 + k + another;})
    .build();

• mesh::sparse_field_builder(obj) - create a sparse field. The syntax is also gridtools::storage-like. The created field shares dimensions with the neighbor connectivity table:

  mesh::sparse_field_builder(obj)
    .type<double>()
    .from<vertex>()
    .to<edge>()
    .initializer([&calc_sign](auto v, auto edge_no) { return calc_sign(v, edge_no); })
    .build();

• mesh::for_each_neighbor(obj, fun, connectivity_or_sparse_field_ptrs...) - iterate over neighbors.

  mesh::for_each_neighbor(
    obj,
    [&](auto index, auto sign) {
       acc += sid::shifted(edge_ptr, edge_stride, index) * sign;
    },
    at_key<connectivity>(ptrs),
    at_key<sign_tag>(ptrs));

Mesh Generator Inputs

• Logical definition of the mesh available from Python. It is enough to export from C++ smth. like:

  vector<array<int, 2>> edge_to_cell();
  vector<array<int, 2>> edge_to_vertex();

• The list of the calls from Mesh API for generated code that where used in the generated calculation(s). Like:

  [conncetivity from vertices to edges]
  [conncetivity from edges to vertices]
  [vertices size]
  [edges size]
  [create edge tmp field with type double and extra dim K of runtime size]
  [create sparse field vertices->edge with type double]
  etc.

  How generated mesh could look like

• for Naive:

  struct my_naive_mesh {};
  
  integral_constant<int, 88> mesh_get_location_size(my_naive_mesh, vertex) { return {}; }
  integral_constant<int, 273> mesh_get_location_size(my_naive_mesh, edge) { return {}; }
  
  auto mesh_get_neighbor_table(my_naive_mesh, edge, vertex) {
    static const int [273][2] table = {{0, 3}, {8, 4}, ...};
    return helpers::neighbor_table_wrapper<edge, vertex>(table);
  }
  // etc.