What is a robust, iterative solver for large 3-d linear-elastic problems?

I'm diving into the fascinating world of finite element analysis and would like to solve a large thermo-mechanical problem (only thermal $\rightarrow$ mechanical, no feedback).

For the mechanical problem, I already grasped from Geoff's answer, that I'll need to use iterative solver due to the size of my mesh. I further read in Matt's reply, that the choice of the correct iterative algorithm is a daunting task.

I'm asking here if there exist any experience on large 3-d linear-elastic problems that would help me narrow down my search for the best performance? In my case, it's a structure with thin, patterned films, and irregularly placed materials (both high-CTE and low-CTE). There are no large deformations in this thermo-mechanical analysis. I can use my university's HPC [1.314 nodes, with 2 AMD Opteron processors (each 2.2 GHz/ 8 cores)].

I think PETSc could contain something interesting, especially the algorithms which do some sort of domain decomposition (FETI, multigrid) but I'm a bit overwhelmed by the options and have no experience. I also like the phrase "geometrically informed preconditioners", but am unsure if this helps me. I have not yet found something focussing on linear continuum mechanics.

Strong scaling (Amdahl) is very important in my application because my industrial partners can't wait a long time for simulation results. I definitely not only appreciate answers, but also recommendations for further reading in the comments.

• Are you solving static problems? Seems yes. For dynamic or time-harmonic problems, the answer would be different, I think. – Hui Zhang May 28 '12 at 15:16
• static yes. Dynamic is too expensive. – Sebastian May 28 '12 at 19:05

Assuming that your structures are actually 3D (rather than only thin features, perhaps discretized with shell elements) and that the model is larger than a few hundred thousand dofs, direct solvers become impractical, especially if you only need to solve each problem once. Additionally, unless the structure is always "close" to a Dirichlet boundary, you will need a multilevel method to be efficient. The community is divided between "multigrid" and "multilevel domain decomposition". For a comparison of the mathematics, see my reply to: What is the advantage of multigrid over domain decomposition preconditioners, and vice versa?

The multigrid community has generally been more successful at producing general purpose software. For elasticity, I suggest using smoothed aggregation which requires an approximate "near null space". In PETSc, this is done by choosing PCGAMG or PCML (configure with --download-ml) and calling MatSetNearNullSpace() to provide the rigid body modes.

Domain decomposition methods offer an opportunity to coarsen faster than smoothed aggregation, thus possibly being more latency tolerant, but the "sweet spot" in terms of performance tends to be narrower than smoothed aggregation. Unless you want to do research on domain decomposition methods, I suggest just using smoothed aggregation, and perhaps try a domain decomposition method when the software becomes better.

• thanks a lot for this very informative answer! What exactly do you mean by close to a Dirichlet boundary? Close in terms of element count? – Sebastian May 29 '12 at 8:32
• Close in terms of distance, measured in elements or subdomains (for one-level domain decomposition with strong subdomain solves), along a path of strong material. If information has to travel through many subdomains to determine a local solution, one level methods will converge slowly. Note that one strong connection is not enough for elasticity, all rigid body modes have to be controlled. – Jed Brown May 29 '12 at 12:09

I would say the canonical choice for this problem would be the Conjugate Gradient solver plus an algebraic multigrid preconditioner. For PETSc, hypre/boomeramg or ML would be the obvious preconditioner choices.

All of these components when used through PETSc scale very well to thousands or tens of thousands of processors if the problem is large enough (at least ~100,000 degrees of freedom per MPI process).

• Note that BoomerAMG (and classical AMG in general) does not use null space information or otherwise ensure that coarse spaces can accurately represent rotational modes. You can (and should) try it, as well as splitting the components and solving them separately (use PCFIELDSPLIT in PETSc), but smoothed aggregation is usually more robust for elasticity. – Jed Brown May 28 '12 at 19:22

You can easily create a rather good preconditioner for elasticity if your material are not too incompressible (Poisson ration $< 0.4$) by using separate displacement components preconditioner. The idea is simply to have a block diagonal matrix with 3 blocks, each of which represents the couplings along the same physical dimension $(x,y,z)$ ie: $K_{xx}$, $K_{yy}$, $K_{zz}$.

In that case you can use less advanced AMG methods to compute an approximation of each block inverse and get a pretty good preconditioner.

Walter Landry developed a code for three-dimensional elastostatic deformation using adaptive multigrid. You can find the code at

https://bitbucket.org/wlandry/gamra

You can include the effect of the thermal forcing with isotropic eigenstrain and equivalent body forces. Once these are in place, the solver would work just fine.