(Proposal WIP) Separate Geometry and Discretization (Enabling of high order finite elements, multi-physics on multi-level)

[rrsettgast]

WIP. This initial post is a rough proposal, and will evolve based on follow-on discussion.

Currently we have a situation where we have mesh levels, which contain the mesh geometry and discretization data. The finite element/volume "discretization" is implicitly assumed to be embedded in the singular geometric "mesh".

           /--- NodeManager
          /      - coordinate data
         /       - field data
MeshLevel --- ElementManager
         \     - quadrature data
          \    - connectivity 
           \--- EdgeManager
            \-- FaceManager

However, in order to effectively support high order meshes, one approach would be to break the MeshLevel up into a geometric BaseMesh consisting of "Cell/Vertex/Face/Edge" data, and a DiscretizationLevel consisting of "ElementManager/NodeManager/FaceManager/EdgeManager". The MeshGeometry object would be pretty static, while the DiscretizationLevel would depend on the order/resolution. The physics solvers would then operate on the DiscretizationLevel (they currently act on MeshLevel), and not the geometric mesh.

To illustrate the relationship between the objects in the proposed BaseMesh and DiscretizationLevel consider the following:

           /--- VertexManager
          /      - coordinate data
 BaseMesh --- CellManager
          \    - connectivity 
           \--- EdgeManager
            \-- FaceManager

Whether or not the BaseMesh object will need to contain sub-Groups such is debatable. I would argue that we could just have the BaseMesh object contain the coordinates of the vertices, and the connectives for the higher topology objects (edge, face, cell), thus flattening the structure.

          /--- NodeManager
         /       - field data
DiscretizationLevel --- ElementManager
         \               - quadrature data
          \              - connectivity 
           \--- EdgeManager
            \-- FaceManager

The NodeManager contains all the "nodal fields", but what to do with the coordinate data that is held in the BaseMesh::Vertex. It makes sense that the Node would have an interface for the coordinates, and it can be a dynamic relationship.

• If the DiscretizationLevel object/refinement is first order (i.e. nodes overlay on vertices) then the Vertex coordinates should be referred to instead of copied.
• If the DiscretizationLevel object is different order/refinement ( #nodes != #vertices ) then the coordinate data may be copied if it makes sense....otherwise a more complex relationship. Summary of options:
  1. NodeManager contains coordinate data. It is distinct from the Vertex coordinate data.
  2. NodeManager doesn't contain coordinates, which means that it is just a container for field data. Access to the BaseMesh ::Vertex is required to calculate the coordinates of a node. This means that outside the context of an element, a node may not be able to know its coordinate. We would want to do this to avoid carrying around a lot of high order/resolution nodal data which may be calculated from much less data dense BaseMesh::Vertex.
  3. We allow NodeManager to contain coordinate data, but only use that data sparingly. We do not refer to NodeManager::coordinate inside of element kernels, instead relying on BaseMesh::Vertex in that context.
• The Cell holds the connectivity with the vertices.
• The Element holds the connective with the nodes.
• The Element holds quadrature data.
• Face/Edge managers hold quantities/field data that are "held" there.

Usage

The only major usage change may be when the current NodeManager requires coordinate data...if we remove coordinate data from NodeManager. Kernels will have be modified to take in BaseMesh::Vertex coordinate data and calculate the coordinates of the nodes. Of course, this can be embedded in the interface for NodeManager, but we will end up with some usage issues to properly move the data between host/device.

Interaction between DiscretizationLevel

DiscretizationLevel will be discretization on top of the geometric "BaseMesh".

• The DiscretizationLevel will have a distinct resolution and/or order from the other DiscretizationLevels.
• A key point is that we must define a restriction/prolongation operation between objects at different DiscretizationLevel s s.t. physics equations can be applied to one DiscretizationLevel , and coupled with physics/fields from another DiscretizationLevel .
  • we may want to have some basic RP operators, but allow a physics solver to override the method to suit the enforcement algorithm.

One example where this will be used the plan for fracture modeling. The Embedded approach has yet to explore the use of a fracture criteria for propagation. The current plan is to have a coarse Embedded representation, and then have a fine scale continuum PhaseField representation, and tie the two Level together. Need Reference.

Communications

This proposal will result in modifications to the ghosting process. Ghosting will have to be setup on the BaseMesh, and subsequently each DiscretizationLevel.

[corbett5]

For first order elements are we going to duplicate the geometric objects?

[rrsettgast]

@corbett5 See changes to top level proposal.

[corbett5]

I would suggest that BaseMesh just be a specialization for a first order DiscretizationLevel.

[corbett5]

There's also been interest in adding support for structured meshes and simplifying the logic to get to a single element and I bring them up here because the solutions all influence each other.

To simplify the element indexing (and I think life in general) I would squash all the element regions and sub-regions into one, storing all the elements in a single allocation (so a single index) and then everything would operate on a set of elements. This logic could be re-used to allow allocating a field on a set of other topological objects as well. If we allow non-contiguous element local indices we can have performant element insertion (for surface generator) by giving each sub-region some extra capacity to grow. Basically an ArrayOf2DArrays where each inner two dimensional array is the element-to-node map for a sub-region.

I think that would pair nicely with a structured mesh implementation as well since if each rank has a ijk structured mesh the complexity of region and sub region is not needed, although operating on a set of elements would most likely be required. I also second @rrsettgast proposal to have a single BaseMesh object instead of separate managers for each topological object since that would be the only sensible way to do it for a structured mesh. As for the complexity of launching kernels on different mesh types we could create a MeshView object that has read only access to only the topological mesh connectivity (not the nodal positions or fields). This object would implement the various topological requests such as elements to nodes, edge to face etc. Since this information should remain unmodified except when using the surface generator it would not incur excess data movement.

[corbett5]

As far as I can tell there are a few use cases

• A n-dimensional field on all elements. This can be done using a single flat n+1 dimensional array like we do with the nodes or if we want to have extra capacity for the "sub-regions" and can live with non contiguous local indices an array of n-dimensional arrays where the outer array index is the sub-region and the inner index gets you the the sub-region local element index. In both ways you could access this data via a single index.
• A n-dimensional field on all the elements where the size of the field varies with the element type. This has to be stored as an array of n-dimensional arrays. For example this would be how the element-to-node map is stored. You could also access this via a single index.
• A n-dimensional field on a set of elements. Here we would have an array1d< localIndex > setIndices that holds the indices of the elements in the set, then we would also have an array< T, N > field that holds the values of the field such that the value at field[ i ] corresponds to element setIndices[ i ]. You would need two indices to access this data (one for the set and one for the element in the set).

As far as I can if we made each region/subregion a set of elements and then allocated everything on a set this would be equivalent to our current interface.

Again I'm not saying we should do this, but I think it's worthwhile thinking about what we would want to do if we were designing it from the ground up and then work backwards from there to adapt it to reality.

[rrsettgast]

To some degree, the ElementViewAccessor and MaterialViewAccessor illustrate the difficulty. We are currently allocating data on the ElementSubRegion, then constructing these accessors s.t. we can access the data above the ElementSubRegion Group through uses of an index for the region and subregion. I think this interface would remain. However, the important thing is that the data could be more naturally accessed outside the subregion context, without having to construct the accessor, as they are implicitly defined. I think shifting the storage location to a legitimate array of arrayofarrays at the ElementManager level would be great, but also retaining the ability of the ElementSubRegion to have a view to the sub-Array (sub-arrayofarrays).

I think we had talked about this a few times a couple of years ago.

[rrsettgast]

To be concise, for n-elements total in, r-regions with nr-elements in each region, rs-subregions with nsr elements in each subregion, we have fields that are:

1. Data allocated for every element. So that is n total entries. Things like volume come to mind. As the first bulletpoint in @corbett5 post implies, these are easily accessible from every level of the element group/set hierarchy using offets and we can include padding allocation for easier insertion.

elementRegionManager->registerWrapper< ArrayXd<real64> > ( "volume" );

Note ArrayXd would have to allow for the gaps if we want them....but I think we would want contiguous access from the top level...so no gaps would be my call. Then when we get it in the scope of the ElementSubRegion we would have to create an upward looking interface... for example,

OffsetArrayView<real64> view = elemSubRegion->getElementReference< ArrayXd<real64> >( "volume" )

This view would have length nsr, while if we got it at the ElementRegion scope, it would have a length of nr.
--> 1. ElementRegionManager sees this as an Array<> of length n,
--> 2. ElementRegion sees this data as OffsetArray<> of length nr,
--> 3. ElementSubRegion sees this data as OffsetArray<> of length nsr.

2. Data allocated for every element in a ElementRegion. Allocations will be on some regions, but not others. Things like the average 'fluidPressure' in the element would be sized in nr length blocks (conceptually)...or constitutive parameter data.

   1. ElementRegionManager sees this as an ArrayOfArraysof size (num region, nr). You still need the region number to get to the data.
   2. ElementRegion sees this data as Array<> of size nr. You don't need the subregion index to get to the data.
   3. ElementSubRegion sees this data as OffsetArray<> of size nsr

3. Data allocated for every element, but varying in dimension across the subregions. The element to node map is the use case here.

   1. ElementRegionManager sees this as an Array1d<Array1d<Array2d<>>>. We would want a new contiguous data type that was contiguous, it would be Array1dOfArray1dOfArray2d<> or something nuts like that. You would still need the region/subregion number to get to the data. Without the fallout of collapsing everything, I don't see a way out of this.
   2. ElementRegion sees this data as ArrayOfArray2d<>. You do need the subregion index to get to the data.
   3. ElementSubRegion sees this data as Array2d<> of size nsr.

4. Data allocated for every element in a subregion, with varying dimension across subregions....like the preceding case...but the data could be higher dimension that 2d (num elem x num quadrature). So this clearly constitutive quadrature data. This is case 3, but worse because of the generalization of the data at the subregion level.

[corbett5]

If we want to be able to use a single index to access volume then it needs to be stored in one single allocation. If we chose to store volume in a standard LvArray::Array< real64, 1 > the adding a new element to a sub-region is a mess, and not just performance wise. Instead of appending you now have to insert. Worse than this though you now have to re-number the local indices of every element in the entire mesh (well you don't have to renumber the elements with lower local indices, but same thing).

If instead we relax things a bit to allow for non-contiguous element local indices then we can store volume in something that is basically just an ArrayOfArrays< real64 > where the first index would just be a combination of the region and sub-region and the second is the sub-region element index. On the region level this would be something that behaved just like an ArrayOfArraysView< real64 >. The value of volume corresponding to the element k would be at volume.offsets()[ k ]. This way creating a new element involves only appending and there is no need to renumber the element local indices. This would make writing a loop over every element a little bit harder, but it is still miles better than the current situation.

Right now loops over all the elements look something like this

arrayView1d< localIndex const > const er = ...;
arrayView1d< localIndex const > const es r= ...;
arrayView1d< localIndex const > const k = ...;
array1d< array1d< arrayView1d< real64 > > > const volumeWrapper = ...;
for( localIndex i = 0; i < totalElements; ++i )
{
  real64 & volume = volumeWrapper[ er[ i ] ][ esr[ i ] ][ k[ i ] ];
}

With non-contiguous local indices you could write the same loop as follows

ArrayOfArraysView< real64 > const volumeWrapper = ...;
auto const volumeView = volume.viewAs1D();
for( localIndex elemIdx : volumeView.validIndices() )
{
  real64 & volume = volumeView[ elemIdx ];
}

When running in parallel either on the host or device it becomes a bit more complicated but still pretty straightforward.

• We could very easily create an abstraction to launch one kernel per contiguous element block. This would at worst launch one kernel per sub-region.
• We could launch a single kernel and then each iteration perform a search on the offsets to map the iteration index to the local element index. My guess is that for non-trivial kernels the overhead here would be minimal.
• We could maintain one single array that maps the iteration index to the element index. Creating a new element could be done by appending a new value and incrementing the entries for all elements with greater local indices than the new element.

Case three seems to be the same as the first case except you instead of using an ArrayOfArrays to store volume you would need an ArrayOf2DArrays to store the element to node map (or obviously the more general ArrayOfNDArrays). You would still need the region and sub-region indices to access the data.

One thing worth noting is that with my suggestion when allocating permuted multi-dimensional data on every element as in case 1 such as cellCenter or rotationMatrix it would not be possible to access this data with just the element index. You would need the full triple.

[rrsettgast]

@corbett5
This sub-discussion is probably worthy of its own thread focused on the element hierarchy and associated data storage....but to continue here for now...I just don't think we care that much about inserting new elements as to motivate adding/keeping the ability to insert an element without being forced to shift. Granted, it is a nice thing to have, but I am not sure we need it. I am all for a single allocation at the top level ElementRegionManager.

You are correct about case 3 being similar to case 1. However the varying geometry on the subregion level (assuming we do not conceptually flatten the entire structure) the data could be represented by an ArrayOfArray<localIndex>. First index is the element index, second is the "element local" node indices. Obviously the values are the nodal indices themselves. This would work but discards the fact that all elements in the subregion have the same...hence your suggestions of ArrayOf2DArray...but then the first index is not the element index. We are back to a double index at this point. One idea would be to store this data in an ArrayOfArray<localIndex> and recognize that the data is contiguous...and overlay a ArrayView2d<localIndex> that was an offset of the ArrayOfArray if we wanted to treat it like we do now for a SubRegion set....but this would remove our ability to permute.

[TotoGaz]

I do not understand a lot about all of this right now, but here is how I would consider the question.

The user writes an xml input, mostly gets VTK/SILO files in return.
In the middle we have solvers and kernels; everything is directed toward providing these solvers with the proper information.
My questions would then be:

• Which exact information do those solvers need to do their duty?
• Which exact information do they need to send back to other simulations or to output files?
• Which exact information do we need to access in order to make the multilevel work?
• ...

If we have a clear view of the flows, it will make drawing connections with other parts of the code easier.

E.g. for FEM, I assume we need a description of the test/basis function, the geometric support (cell). Based on this we'll be able to assemble the matrix. Do we need something else and where? Also how do we send the data back to the user?…
Knowing this will help us defining the contracts of the classes we want to create, eventually we'll find out how to implement them.

One convenient way could be to write pseudo code so we know the algorithm and what it needs from which or which element.