Which preconditioning for large linear elasticity problem?

Question

The problem I want to solve is the displacement formulation of the linear elasticity : $$ \nabla \cdot \sigma = 0 \quad \text{in} \quad \Omega \\ \sigma = \lambda ( \nabla \cdot u ) I + \mu (\nabla \cdot u + \nabla \cdot u^t) \quad \text{in} \quad \Omega \\ u = 0 \quad \text{over} \quad \Gamma_D \\ \sigma \cdot n = g \quad \text{ over } \quad \Gamma_N $$ with $u$ the displacement vector, $\sigma$ the stress tensor and $\epsilon$ the deformation tensor.

The weak formulation is the standard one : $$ a(u,v) = \int_{\Omega} \lambda \; \nabla \cdot u \nabla \cdot v + \int_{\Omega} 2 \mu \; \epsilon (u) : \epsilon (v) \\ f(v) = \int_{\Gamma_N} g \cdot v $$

Then there is the FEM discretization over a tetrahedral or hexahedral mesh using the classic lagrange functions basis of degree 1 or 2. (P1 or P2 for tet mesh and Q1 or Q2 for hex mesh)

I can assemble the linear system $ A x = B $ without difficulties but I am unable to solve it for non trivial problems.

I have reduce my test case to the following :

$$ \Omega \; \text{ is the unit cube } \\ u = 0 \quad \text{ for } z = 0 \\ \sigma \cdot n = -0.5 \sin(\pi x) \sin(\pi y) \; e_z \quad \text{ for } z=1 $$ The mesh is a regular subdivision in tetrahedra (or hexahedra) of the cube.

I am using the Conjugate Gradient to solve the linear system. I tested with Jacobi or Gauss-Seidel preconditioners but the CG solver is faster without them.

For Q1 function basis :

• The CG solver does not converge (almost constant residual) when there are more than 300k unknowns (100k vertices)

For Q2 basis :

• It starts to be slow at 200k unknowns and after it won't converge at all

Using Jacobi or Gauss-Seidel preconditioner does not improve the convergence.

For a comparison view, I tested with the heat equation (scalar field) and a similar problem (unit cube, mixed boundary conditions), the CG solver converge quickly with 600k unknowns. As the difference is mainly scalar vs vector field, I guess I am missing something in the preconditioning of vector problem.

A linear elasticity problem with 200k unknowns (8000 hexes, 9000 vertices, Q2 basis) and a trivial geometry seems pretty small but I am failing at it. And the simpler preconditioners do not improve things. So my questions are :

• How should I choose my preconditioners for linear elasticity ?
• How people solving large elasticity problems (millions of unknowns) deal with this kind of issues ?

I would be very grateful if you have any advices or references on the subject, thank you.

Answer 1

You can always try AMG as one preconditioner.

There are several other methods that are used. Look into papers by Johanes Kraus, Neytcheva, Owe Axelsson.

If you use separate displacement ordering you can try Schur complement type preconditioners which are very effective.

You may notice that for Incompressible material $\nu = 0.5$ you get $\lambda \rightarrow \infty$ which is a problem.

What people do is they introduce another equation like this $$ \mu \nabla \cdot u - \frac{\mu^2}{\lambda} p = 0 $$

You use stable pair Q2-Q1 for $u$ and $p$. From this you get a two by two matrix $A$ which you can solve with Schur complement.

I have been working on this for some times now, look for papers by Ali Dorostkar (which is I) and Maya Neytcheva.

Update

There is one thing I forgot to mention here. From Korn's inequality, we know that the block diagonal part of matrix $A$ is spectrally equivalent to $A$ so you can use that as your preconditioner. Meaning that you create a preconditioner $$ D = \begin{bmatrix} A_{11} \\ & A_{22} \end{bmatrix} $$

Since $Ax$ is a vector (lets name it $v$) in each iteration you apply your preconditioner you have $$ D^{-1}Ax = D^{-1}v $$ but you don't want to invert $D$ so what you do is that you solve another system like this

$$ Dy = v $$

now since D is block diagonal you just need to solve two smaller systems. Each of those smaller systems you precondition with AMG.

I know this is complicated but I thought maybe this interest you.