7
$\begingroup$

Ordering of elements in an unstructured mesh is undoubtedly very important for the performance of computations. For example, it determines the structure of sparse matrices arising from PDE discretizations, which affects the performance of most linear algebra operations like matrix-vector product.

Recently we used an unstructured tetrahedral mesh generated by COMSOL Multiphysics. For given mesh, we construct a matrix with the following structure: each row corresponds to some facet $F$ of the mesh and the non-zero elements correspond to the facets of the one or two cells adjacent to the facet $F$, i.e. for tetrahedral mesh each row has at most 7 non-zero elements. The matrix structure for the mesh generated by COMSOL is shown on the following figure.

enter image description here

Simultaneously, we experiment with other mesh generators, but we currently don't have a good way to order the mesh. Using a simple approach of an in-order traversal of an octree of the cell centers, we get the following structure of the same matrix as above, which is not as "nice" as the original:

enter image description here

Unfortunately the COMSOL documentation does not say anything about mesh ordering algorithms. Does anybody have an idea which algorithm it might be using or which algorithms might generate similar (or better) structure?

Improve this question
$\endgroup$
3
  • 7
    $\begingroup$ Just eyeballing the result, that looks like the reverse Cuthill McKee ordering. I suggest, however, that you look at the Deal.II page on mesh orderings before you choose the one to use. The orderings attempt to achieve different goals (bandwidth reduction, Cholesky fill-in reduction, ...) that are not simultaneously achieveable. dealii.org/8.5.0/doxygen/deal.II/namespaceDoFRenumbering.html en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm $\endgroup$
    – Tyler Olsen
    Commented Apr 30, 2017 at 20:25
  • $\begingroup$ Thank you, that indeed looks like the Cuthill-McKee ordering according to the figures at Wikipedia. I'll try to hook the boost implementation of the method into our code and look at the results, although comparing with the Deal.II pictures, the result will most likely depend on the way the algorithm orders vertices with the same degree. Out of curiosity, are there any heuristics for that? $\endgroup$
    – Jakub Klinkovský
    Commented Apr 30, 2017 at 20:55
  • $\begingroup$ RCM is a variant on a breadth-first search, so you really only need to define a tie-breaker for equal-degree neighbors. Not being an expert at all on the nuances of these choices, my first thought would be to first try nothing (i.e. order in which you see them is order that they go in). You might also try something like computing the sum of the node degree + neighbor degrees, which will hopefully break a few ties. Fortunately, these sorts of things are pretty easy to try out with a baseline implementation in place. $\endgroup$
    – Tyler Olsen
    Commented Apr 30, 2017 at 21:02

2 Answers 2

Reset to default
8
$\begingroup$

The idea of "ordering the nodes" in a finite element mesh to improve the computational time of the sparse solver originated in the large structural analysis FE codes of the 70's. Those codes typically used banded or variable-band storage schemes for the sparse matrices so reducing the bandwidth was the main criterion. That is the origin of the old Cuthill-McKee and Reverse-Cuthill-McKee bandwidth reduction algorithms.

Those concepts are essentially obsolete today in all FE codes, such as COMSOL, that use modern sparse solvers.

You say that your second ordering is "not as nice" as your first, presumably because the bandwidth is larger than the first. Bandwidth is a completely irrelevant metric in all modern sparse solvers.

Here is a specific example that you may find surprising. In a classic paper from 1973, nested dissection, Alan George showed that the optimum numbering for the Laplace equation discretized with a simple square finite element mesh was obtained by what he called "nested dissection." The non-zeros for a $7\times7$ mesh with equations ordered according to the nested dissection algorithm are shown in the figure below. enter image description here

As you can see, the bandwidth of the matrix is nearly equal to the number of equations. This is the equation ordering that minimizes the number of floating point operations during factorization.

All modern sparse solvers that I am aware of include algorithms for reordering the equations for computational efficiency and they are invoked by default. Many include multiple algorithms and will try two or more to find the one which minimizes the factorization time.

The latest COMSOL documentation I have found shows that the direct sparse solver options are MUMPS, PARDISO, and SPOOLES. All three of these have excellent algorithms for reordering equations. They include algorithms that attempt to find an ordering similar to George's nested dissection on general meshes.

For your specific problem, it is worthwhile to review COMSOL's guidelines regarding choice of direct solver. And you might want to try one of the iterative solver options where the equation ordering is irrelevant. But, in either case, the numbers you use to label the nodes in your mesh are not going to have any significant effect on the computational time. So pick ones that are most convenient for you.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Obviously the optimality of the ordering depends on the method that is used to solve the system. For direct solvers the ordering usually involves minimizing fill-in, which is irrelevant for iterative solvers. But I wouldn't say that ordering is completely irrelevant for iterative solvers. For example, the dominant operation of Krylov methods is usually the matrix-vector multiplication. If the matrix is so large that the vector does not fit into cache, then the cache misses due to bad ordering will kill the performance of I'd say any Krylov method, even GMRES where SpMV is not dominant. $\endgroup$
    – Jakub Klinkovský
    Commented May 1, 2017 at 7:40
1
$\begingroup$

Based on the comment by @TylerOlsen, I've applied the boost implementation of the reverse Cuthill-McKee ordering on the same matrix from the original question and here is the result:

enter image description here

The structure looks very similar, but the bandwidth is significantly smaller compared to the original matrix with COMSOL ordering. Also changing the starting vertex in the Cuthill-McKee algorithm (e.g. enforcing to use the same starting vertex as COMSOL) does not have a significant effect on the bandwidth for this particular mesh. I consider this to be even better result.

Since I can't formally accept comments, I'm going to self-accept this answer, but of course the credit goes to @TylerOlsen.

Improve this answer
$\endgroup$

Your Answer

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

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