I am working on modeling incompressible elasticity at finite strains.

$$ \mathrm{Div} \boldsymbol P = \boldsymbol 0, \quad \boldsymbol X \in \Omega_0 \subset \mathbb R^3, \\ J = 1, \qquad \boldsymbol X \in \Omega_0 \subset \mathbb R^3, $$

with given Dirichlet BCs

$$ \boldsymbol u = \bar{\boldsymbol u}, \quad \boldsymbol x \in \partial \Omega_0. $$

Here, $\boldsymbol P(\boldsymbol F)$ is the Piola stress tensor, $\Omega_0$ the material configuration, and $ J = \det \boldsymbol F $ is the determinant of the deformation gradient.

The incompressibility constrained $ J=1$ is incorporated via the Lagrange-Multiplier method leading to a saddle-point problem $$ \begin{bmatrix} \mathbf A & \mathbf B^T \\ \mathbf B & \mathbf 0 \end{bmatrix} \begin{bmatrix} \boldsymbol u \\ \boldsymbol p \end{bmatrix} = - \begin{bmatrix} \boldsymbol r_u \\ \boldsymbol r_p \end{bmatrix}, $$

which must be solved at every Newton iteration. So it may happen that this linear system must be solved 50-100 times to solve the above BVP, depending on the applied load.

$\boldsymbol A$ is symmetric and positive definite and $B_{ij}^t = \int_{\Omega_0} \left[\boldsymbol F^t \cdot \boldsymbol v_i \right] : \left[ \boldsymbol C^{-1} J q_j \right] dV $ with $ \boldsymbol C = \boldsymbol F^t \cdot \boldsymbol F$, that is, also symmetric and positiv definite. $\boldsymbol v$ and $q$ are the test function associated with the displacment and the pressure (Lagrange multiplier), respectively.

Note that $ J = 1$ is equivalent to $ \mathrm{div} \ \boldsymbol u = 0$ in small-strain elasticity, in which case the PDE is quite similar to the Stokes equation. There are a lot of conversations here talking about preconditioners for the Stokes equation, but I have not found discussusions for the non-linear regime.

In 2d for moderate problem sizes (less than 100k degrees of freedom), I am aware that a direct solver like UMFPACK is the way to go. In 3d, I expect problem sizes less than 500k DoFs. Clearly, a direct solver is no longer applicable and I have to come up with an iterative solver along with a preconditioner.

How could a good preconditioner look like? Can the techniques applied to the Stokes equation, such as the Schur complement, be used here too?


2 Answers 2


There is no difference between the linear and nonlinear Stokes problems as far as preconditioning is concerned. That is because at the end of the day, you always have to linearize the nonlinear problem (e.g., via a Newton iteration) and then solve the linear problem -- which is conceptually a Stokes equation, though perhaps with spatially variable coefficients. The right preconditioner for these linearized steps are then simply the preconditioners that have already been developed for either the Stokes equations, or for Stokes equations with variable viscosity.

  • $\begingroup$ So you would say that my B-matrix has the structure of a Stokes equation? If so, is the "Block Schur complement preconditioner" you describe in step-22 a good approach for my problem? In this question, scicomp.stackexchange.com/questions/33859/…, the block-proconditioner looks different, for instance. $\endgroup$
    – Simon
    Dec 12, 2023 at 8:38
  • $\begingroup$ This question may be specific to deal.ii but fits well in this discussion. I run step-56 with Q1-P0 (constant pressure) element. At refinement cycle 1 (4096 cells, dofs = 14739+4096), it took 311 FGMRES iterations and 2899+322 iterations for approximating A and S inverse, respectively. Using the Taylor-Hood Q2-Q1, there are in total only 223 iterations. Why is the block-preconditioning implemented there not appropriate for Q1-P0 element? I do not want to use quadratic elements for displacements. $\endgroup$
    – Simon
    Dec 12, 2023 at 16:23
  • $\begingroup$ I didn't write step-56, so I can't say why they chose the preconditioner they use :-) But in general, I always recommend the Sylvester-Wathen preconditioner for Stokes-like equations, see the discussion at the end of step-22. $\endgroup$ Dec 13, 2023 at 1:59
  • $\begingroup$ In step-56, they implement exactly the Sylvester-Wathen approach you describe at the end of step-22 (+ the possibility to use AMG as preconditioner). So the issue is clearly related that one uses Q1/P0 instead of Q2/Q1. Shall I open a question on the list or is it clear what I mean? $\endgroup$
    – Simon
    Dec 13, 2023 at 5:57
  • $\begingroup$ So my question can be reformulated as "why is the preconditiong at the end of step-22 so inefficient if one uses Q1/P0 instead of Q2/Q1?" $\endgroup$
    – Simon
    Dec 13, 2023 at 5:59

Although the matrix structure looks the same, solving the matrix system for the Stokes problem is not the same as solving the matrix system for incompressible hyperelasticity problems. The condition number of matrices would be different, changing with iterations in the Newton-Raphson loop for the hyperelastic problems.

Nevertheless, we can still use the same approach of block preconditioners and approximated Schur Complement as in the Stokes system, but with different types of Krylov solvers and preconditioners. In my experience, the preconditioners for Stokes and Navier-Stokes do not work the same for incompressible hyperelastic problems because of the severe mesh distortions. It depends on a lot of factors, element aspect ratios, mesh distortions, and material properties.

Check example 19 on MFEM website. This is exactly what you are looking for.

  • $\begingroup$ Out of curiosity, why does the condition number changing between iterations matter? I mean, of course it does, the nonlinear coefficients change. But you of course also update the preconditioner accordingly. $\endgroup$ Dec 24, 2023 at 2:53
  • $\begingroup$ @WolfgangBangerth, it is not the change of condition number that matters most, but the differences in condition number itself between the Stokes system and system for incompressible hyperelasticity. $\endgroup$
    – Chenna K
    Dec 28, 2023 at 0:20
  • $\begingroup$ @ChennaK The implemented preconditioner in MFEM looks very similar to the preconditioner applied to the Stokes system (see e.g the linked step-56 tuturial above). Minor differences are the outer solver (GMRES instead of FGMRES), the approximation of the stiffness block (CG vs GMRES), and a constant weighting factor on the mass matrix. Other than that, I do not see any differences. $\endgroup$
    – Simon
    Jan 9 at 13:27
  • $\begingroup$ @Simon, yes. It looks very similar to the one for the Stokes system. You will notice the differences once you run the code. The default settings do not work for refinement level 2 or higher. You can, however, experiment with different solvers, preconditioners and the parameters in the context of hyperelasticity. $\endgroup$
    – Chenna K
    Feb 7 at 10:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.