# implicit vs. explicit domain decomposition methods

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?

• I think for stencil-wise (e.g. FD) assembly, nodal partition is natural. For element-wise assembly, element partition is more natural. But in the actual implementation, one always needs a disjoint partition of d.o.f to keep each d.o.f uniquely valued. To provide a local (on this processor) element partition, in PETSC, you can use AO to set up LocalToGlobalMapping for Mat. Another way is to use DG method, then no d.o.f is shared by two elements. Feb 14, 2013 at 18:07

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.

• Could you elaborate on why an element-wise partition yields worse load-balancing and more communication? With a node-wise partition, you have to send values at the ghost nodes, and with an element-wise partition you have to coordinate values at the shared nodes; I don't see how there's a big difference in communication volume between the two approaches. Nonetheless, your point about greater facility with overlapping Schwarz and AMG is well taken -- I hadn't considered those methods. Feb 14, 2013 at 18:52
• Consider a corner point that is split between $k$ subdomains. Since you have to make the vector consistent before the next iteration, you somehow have to do all-to-all communication between those $k$ parts. If you had a vertex partition, that point would be owned by exactly one process and communication is simpler. It's not a big effect. At the end of the day, the linear algebra only needs vertices. If your computation uses both, then you have to load balance both (e.g., hypergraph partitioning) which is harder than just partitioning vertices. Feb 14, 2013 at 22:49

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.

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.

• What method has supplanted DD then? (Not being snarky here, I honestly don't know.) Feb 14, 2013 at 19:02
• Methods that treat the PDE as one, global problem. You then get one very large linear system that is stored in parallel by letting each processor "own" a certain part of the rows of the matrix and vectors. You then apply a Krylov space method to it which requires doing matrix-vector products with this globally stored linear system. Feb 16, 2013 at 12:02