13
$\begingroup$

I'm interested in an advice for efficient data structures for cell browsing in unstructured cell-based finite volume CFD.

One example that I encountered (in dolfyn cfd code) goes like this (I'll show relevant segment) \begin{listing} do ip=1,Ncel ... do j=1,NFaces(ip) k = CFace(ip,j) ipp = Face(k)%cell1 inn = Face(k)%cell2 if( inn > 0 )then ! internal \end{listing} So we have an array NFaces where the number of faces for each cell is stored. Then CFace array which maps local-to-cell face number to global face number.

The code is face based therefore there is a face data type which stores the serial number of two cells it lies between Face(k)%cell1 and Face(k)%cell2.

Any comments on this or suggestions for alternative approach are welcome.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
10
$\begingroup$

The structure you show is a common choice and equivalent to storing cell-face adjacencies in a CSR matrix format, with the boundary ghost cells in a special place. However, note that FV methods can also be formulated to consist entirely or almost-entirely of face traversal where each face is only visited once (reconstruct to face centroid/quadrature point from both sides, solve Riemann problem, distribute flux back into residual on cells). You can "fake" that by using your cell-based traversal and skipping any two cells that are below the "diagonal" in the sparse matrix, but a popular alternative is to store (leftCell, rightCell) = support(face), in which case faces become first-class entities. That is useful because you typically need a place to store face quadrature points (centroids), face normals. You can also put reconstruction parts (like least squares) into the face-based data structures. The face traversal is seemingly vectorization-friendly because all the sizes are regular, but on the other hand, there are overlapping outputs so you need to organize the traversal to avoid putting it in an inner loop. With this more face-oriented data structure, it is natural to order the face numbers so that each boundary condition type can be applied using a contiguous traversal of faces (also vectorization-friendly).

If you choose this data structure, remember to sort the faces so that the traversal reuses cell data in cache as much as possible. See any of the PETSc-FUN3D papers for performance analysis of face ordering and related optimizations.

Improve this answer
$\endgroup$
2
  • $\begingroup$ If you loop over the faces, you will need to get information from a leftCell and rightCell to compute the fluxes, say face 1 has leftCell 1 and rightCell 10, face 2 has leftCell 6 and rightCell 31, ... thereby jumping through the memory. How would that be vectorization friendly? $\endgroup$
    – chris
    Nov 29, 2012 at 20:17
  • $\begingroup$ As mentioned above (and discussed in the PETSc-FUN3D papers), you order the faces to reuse cache. The result is like a "one-sided" cell traversal in which each face is visited only once. $\endgroup$
    – Jed Brown
    Nov 30, 2012 at 4:58
3
$\begingroup$

I know this question is already answered, but here's a similar single-face based looping storage which is implemented in the OpenFOAM C++ library:

Each cell has an index (ID) in a cellList. Two lists are defined for all faces: "face internal owner" and "face neighbour". The length of both face lists corresponds to the number of internal faces in the mesh. A face owner will be the cell with the lower ID in the cellList (opposite for face neighbour). Boundary faces are written last, and they have outward oriented normals (from the solution domain), and of course, only one owner cell. The face area normal is oriented so that it looks outwards from the owner cell to the neighbour cell.

This works well for e.g. flux calculation. The flux is evaluated once per face, and it is added to the sum of total faces for the owner cells, and deducted from the neighbour cells (the summation/deduction is decided based on the orientation of the face area normal). The boundary faces are sorted and stored at the bottom of the face list, allowing boundary conditions to be defined as slices of the face list (begining label, end label of the boundary patch),simplifying thus the implementation of the boundary conditions, as well as inhancing efficiency of the updating process for the boundary conditions, since it is relying on the solution provided by the operations on internal faces.

Since the boundary faces are agglomerated into patches, the inter-process communication is defined for coupled (processor) patches, and pre-defined. This means that as soon as there is a loop over the boundary mesh, the top level access functions envoke wrapped MPI calls, making such code "automatically" parallelized, if it relies on the above explained face-based connectivity.

Improve this answer
$\endgroup$
3
  • $\begingroup$ No problem, I'm glad to see that this description is useful to someone.. :) Are you working with OpenFOAM as well? $\endgroup$
    – tmaric
    Dec 6, 2012 at 11:17
  • $\begingroup$ I used to, a little bit in the past. I generally tend to stay away of accepted trends and try to reinvent the wheel. That's my Tao. $\endgroup$
    – Johntra Volta
    Dec 6, 2012 at 11:37
  • 1
    $\begingroup$ Your Tao is the opposite of the Tao of Computer Science: "Don't re-invent the wheel". But I can understand it, it's appealing to do stuff from the scratch! :) $\endgroup$
    – tmaric
    Dec 6, 2012 at 12:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.