exasim-project / NeoFOAM

WIP Prototype of a modern CFD core
20 stars 1 forks source link

Feature adjacency, for mesh connectivity - stencils. #7

Closed bevanwsjones closed 6 months ago

bevanwsjones commented 7 months ago

Disclaimer: Feel free to close this if it's not needed, I am just putting it up 'early' to not spend too much time on it if it's not wanted/the wrong approach (I tried to understand OpenFOAM's approach - but still a bit blurry to me). Also, I started the work on the blas-feature branch, so it has a few commits from that branch - my bad :disappointed: . I have also updated some of the primitives, I wanted to add some syntactic sugar just to make it easier to use, again if this needs to be removed please let me know.

The central feature is the deviceAdjacency, which takes its name from the mathematical concept of an adjacency graph. These graphs make up the backbone of CSR/CSC matrix storage, where the matrix data is stored in three contiguous arrays/vectors, here a Kokkos::Veiw. One for the floating point data, one for the associated floating point's column index, and one for the offsets per row (in the context of a CSR matrix). Unlike a CSR matrix, the deviceAdjacency has no 'floating point' container, and the values stored in the column vector are the indexes indicating the first-order neighbours in the graph.

In an undirected graph (symmetric CSR matrix) the graph represents connections between the same mesh entities, e.g. cell-cell, face-face, etc. By extending the class to allow for undirected configurations, the class can be used to store the connections between different mesh 'entities', cell-face, cell-vertex, vertex-face, etc.

Typically these types of containers are very quick to access, but slow to build compared to the 'vector of vectors' approach to connectivity, however, you are guaranteed good cache hitting for sequential lookups and no memory hops' which may occur in a vector of vectors layout (especially if connectivity changes during run time with no management). Further, the reason for not going for Kokko's multidimensional view is due to them having to be 'square' (fixed cardinality in each dimension), which is not the general case for CFD unstructured meshes. Further, the library discourages an equivalent 'vector of vectors' approach, a 'view of views'.

Example Usage:

NeoFOAM::localAdjacency<false> cell_cell; // undirected graph
cell_cell.insert(1, 2); // adds connection between cell 1 and cell 2.

localIdx x = cell_cell(1)(0); // x = 2    - we can try use the subscript operators[], I prefer the square brackets tbh.
localIdx x = cell_cell(2)(0); // x = 1
size_t cell_cell(0).size(); // = 0
size_t cell_cell(1).size(); // = 1
size_t cell_cell(2).size(); // = 1

NeoFOAM::localAdjacency<true> cell_face; // directed graph read as cell-to-face connectivity.
cell_face.insert(1, 2); // adds a connection between cell 1 and face 2 (importantly not the reverse).

localIdx x = cell_face(1)(0); // x = 2
size_t cell_face(0).size(); // = 0
size_t cell_face(1).size(); // = 1
size_t cell_face(2).size(); // = 0

This still to do:

Finally, I am struggling with cuda enabled Kokkos, more on the Cuda side, so I have not started on benchmarking.

HenningScheufler commented 7 months ago

Thanks Bevan it look really good! I will merge it into proof of concept and will try to test it with the gradientOperator,

The face based approach requires atomics that where about a factor 2 hit on performance and cell_face stencil atomics are not required.

bevanwsjones commented 7 months ago

So I am going to slim down the requirements for this pull request in the interest of speed. I will open an issue for the remaining items - I have added them to the class todo list.

Currently the only remaining issue I think there is, is to ensure the constructors check and correct adjacency ordering. From there I think it's pretty much good to go, unless there are other points. :) Also need some guidance on commenting - now that we using Sphinx etc.

New issue tasks (after this pr):

bevanwsjones commented 6 months ago

One or two further things to sort out:

  1. What is the approach to error handling, do we throw or terminate and do we have a set of 'maco' like functions to aid in standardising the approach?
  2. How do we death test? Currently, I have not added death tests because gtest complains (understandably so).
  3. What is the approach to testing considering test coverage good enough? I have checked the 'testing' checkbox complete, but there definitely could be further coverage.
HenningScheufler commented 6 months ago

can you change the pull request so it merges into proofofconcept?

bevanwsjones commented 6 months ago

can you change the pull request so it merges into proofofconcept?

is this correct - or does it need to go to GPL?