implicit vs. explicit domain decomposition methods

Question

I've been working on a finite element code on unstructured methods, which I've parallelized using the Schur complement method. Here's a summary of how I did it:

1. Assign each triangle of the mesh to a domain
2. For each node, determine which domain it is in or whether it is on the artificial boundary
3. Send each processor its part of the mesh -- the triangles it owns as well as any nodes it owns, even if it shares those nodes with some other domain
4. Build the matrices, in parallel, on each processor
5. Solve some stuff, using some way of exchanging node data between processors

I'm solving a quasilinear system of elliptic PDE, so each linear solve is just an iteration of Picard's method. In particular, I need to change the entries of the stiffness matrix at every step, since it depends on the gradient of the solution. (The non-zero structure of the matrix doesn't change.)

This all works just fine. However, it disagrees with the usual approach I see taken.

My code divides up the triangles between processes, and some of them have to share nodes. This means a extra storage, but the entries of the stiffness matrix can be filled completely in parallel. On the other hand, libraries like PETSc and METIS divide up all the nodes into disjoint sets and send them, along with their corresponding rows in the stiffness matrix, to each process. But, you have to communicate to fill the stiffness matrix.

So: why do the big scientific libraries have a marked preference for dividing up the nodes instead of the elements? Am I missing something?

Answer 1

Partitioning by elements leads to minor simplifications in assembly at the expense of poor load balancing in the solver, inability to use many popular solution methods, and more communication to apply the operator. The most common reason to use an element partition is simply because the code to handle overlap or to communicate contributions into off-process rows during assembly has not been written. (PETSc will do the communication automatically.) It's rarely faster in the long run.

Now if you only assembly locally, you end up with "Neumann problems" on each subdomain. Since these problems have no boundary conditions, they are singular for differential equations with no zeroth order "mass" term. This requires use of much more complicated preconditioners, such as FETI-DP or BDDC, for which the performance is quite sensitive to the partition and that usually require special work if the character of the physical system changes. If you assemble to a vertex-based operator (like PETSc's default matrix formats), you can use many methods that are "easier" to use, such including additive Schwarz and algebraic multigrid. (Broadly speaking, "easier" means that they have a broader performance "sweet spot", do not finnicky null-space issues in which local mesh problems can cause complete breakdown, have better diagnostics, and require less customization for new problems.)

Some people, including the 2004 Gordon Bell Prize winner, use a vertex partition everywhere and integrate redundantly in the overlapping elements to assemble the operator without communication. That is usually best for performance, but requires that the mesh be extended to include one cell of overlap, which may or may not be convenient in your code.

Answer 2

Both are valid techniques with their own tradeoffs. PETSc, in particular, has a history as a linear algebra library, and given that, it makes sense for it to assign rows distinctly to processors. For finite element methods, which have element matrix formation, and often some sort of matrix assembly, it often makes sense to assign elements to uniquely to processors. Some codes do element matrix formation this way, then choose an assignment of nodes to processors before calling PETSc during matrix assembly and solution.

All of these techniques can be made to work.

Answer 3

Part of the disconnect you see is historical preference. Domain decomposition methods were popular in the 1990s when there were no linear algebra libraries around that handled the (rather complex) task of communicating matrix entries if you split vertices (or, rather, degrees of freedom) into disjoint sets; processor numbers on parallel machines were, at that time, also in the few dozens. If you had to write algorithms for parallel communication of data yourself on such small (by today's standards) machines, domain decomposition methods made good sense.

But times have changed. Libraries such as PETSc handle the communication of matrix elements just fine, without user interaction. And machines today have hundreds, thousands, or more processors, and domain decomposition methods simply do not scale well to these kinds of machines. Consequently, domain decomposition methods that split on cells are not all that popular any more.