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:Assign each triangle of the mesh to a domainFor each node, determine which domain it is in or whether it is on the artificial boundarySend 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 domainBuild the matrices, in parallel, on each processorSolve some stuff, using some way of exchanging node data between processorsI'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?
implicit vs. explicit domain decomposition methods
linear algebra;finite element;partitioning;domain decomposition
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.