{'fraction_non_alphanumeric': 0.030227988479448735, 'fraction_numerical': 0.0093133428440688, 'mean_word_length': 4.662551440329218, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parallel inexact Newton-Krylov and quasi-Newton solvers for nonlinear elasticity

In this work, we address the implementation and performance of inexact Newton-Krylov and quasi-Newton algorithms, more specifically the BFGS method, for the solution of the nonlinear elasticity equations, and compare them to a standard Newton-Krylov method. This is done through a systematic analysis of the performance of the solvers with respect to the problem size, the magnitude of the data and the number of processors in both almost incompressible and incompressible mechanics. We consider three test cases: Cook's membrane (static, almost incompressible), a twist test (static, incompressible) and a cardiac model (complex material, time dependent, almost incompressible). Our results suggest that quasi-Newton methods should be preferred for compressible mechanics, whereas inexact Newton-Krylov methods should be preferred for incompressible problems. We show that these claims are also backed up by the convergence analysis of the methods. In any case, all methods present adequate performance, and provide a significant speed-up over the standard Newton-Krylov method, with a CPU time reduction exceeding 50% in the best cases.


Introduction
Nonlinear elasticity is a continuum framework for modeling the deformation of an elastic body, and has a large spectrum of applications, ranging from industrial to academic. The flexibility of continuum mechanics resides in its abstract formulation starting from fundamental principles, which allows for an accurate representation of many different materials through constitutive modeling and complex boundary conditions [Hol02]. After a precise physical phenomenon has been devised and modeled, the resulting equations yield a complex system of nonlinear Partial Differential Equations (PDEs). Limited progress has been achieved in terms of its numerical analysis [CD04, ADVL + 13], so it is difficult to know a-priori efficient strategies for its numerical approximation. Indeed, the problem is not convex but only polyconvex [Cia21], which deteriorates the performance of many well established numerical methods. The gold standard is Newton's method [ZTNZ77], consisting in a second order approximation of the associated variational principle (or a first order approximation of the Euler-Lagrange equations, referred to as the Newton-Raphson method). This method is popular mainly due to its robustness and low iteration count, as it converges quadratically whenever a good initial guess of the solution is considered [WN99]. Its main drawback is that it requires the repeated assembly of the Jacobian matrix, associated to the tangent problem, which usually becomes a bottleneck of the solution process (see [JHN11] for an example in CFD). This can be more evidently appreciated when using higher order approximations, required for avoiding numerical locking effects [EG13].
In general, not much attention has been given in the community to other methods, such as linearly convergent (gradient descent [Rud16], Richardson [Saa03]) and superlinearly convergent (quasi-Newton [WN99], inexact Newton [DES82]) ones. These methods yield a higher iteration count due to their reduced convergence rate, but can potentially present a drastically overall reduced computational complexity. Indeed, the use of different nonlinear solvers has shown enormous speed-ups in other problems [SHOL06,BRKN17], where the properties of the equations have been exploited to devise a more adequate nonlinear solution method. Some recent works have numerically explored quasi-Newton methods for the mechanics [LLS05,GP88,LBK17], with their main focus being on the nonlinear iterations and CPU time. This type of study does not take into account the linear system involved in each iteration of the method, and thus is not conclusive with respect to the applicability of these methods in an HPC setting, where problems with millions of degrees of freedom need to be solved, possibly for many time steps. A numerical method can be deemed adequate in such case if at least is satisfies the following requirements: i) the nonlinear and linear iterations incurred during the solution procedure are independent of the mesh size used, ii) the method is robust with respect to the model parameters, at least in the scenarios of interest and iii) the method is (strongly) scalable, meaning that the overall solution time improves when more processors are used. We finally highlight [WDE07], where a novel conjugate Newton method was proposed and analysed for nonlinear elasticity.
The scope of this work is to provide a first step in the adequate usage of alternatives to Newton's method, where we focus in superlinearly convergent methods as they provide an excellent overall computational complexity [WN99]. We stress that all gradient descent algorithms we tested on a preliminary phase failed, and this is indeed consistent with the literature: there are no works, up to our knowledge, of descent nor fixed point algorithms for nonlinear mechanics. We focus on three tests: i) Cook's membrane test, an almost incompressible problem, ii) a twist problem, where we test an incompressible material and iii) a cardiac modeling test, where we use an idealized left ventricle geometry to model a human heartbeat. The work is structured as follows: in Section 2 we review the mechanics problem and fix some notation, in Section 3 we describe the nonlinear solvers to be used throughout this study, in Section 4 we define the tests in which the methods will be tested together with a thorough description of the linear solver involved, in Section 5 we recall the convergence results of the methods and discuss what their expected performance should be, in Section 6 we show and comment the results obtained and we conclude with a discussion of the results in Section 7.

The hyperelasticity problem
Consider a connected domain Ω ⊂ R d∈{2,3} that represents the geometry to be deformed, and define its Dirichlet and Neumann boundaries as ∂Ω D and ∂Ω N respectively, such that ∂Ω = ∂Ω D ∪ ∂Ω N . We look for a displacement d : [0, T ] × Ω → R d such that it solves the momentum conservation equation where d 0 is the initial displacement, v 0 is the initial velocity,t is the surface traction, d D is a distributed load, () denotes a time derivative, and P is known as the Piola stress tensor, usually obtained from a predefined Helmholtz potential Ψ such that where F = I + ∇d. The existence of such a potential is actually a hypothesis, and whenever it exists we refer to the material as a hyperelastic material. Approximating this problem by neglecting the inertial term d results in the well-known static mechanics, or quasi-static whenever one of the problem's data is timedependent. If we require the volume to be conserved through J := det F = 1, we obtain an incompressible elasticity problem. This can be instead approximated by further penalizing the deviation of J − 1 from 0, which results in a modified potential where the isochoric componentF = J −1/d F is such that detF = 1. This separates the energetic contribution of volumetric deformation, and is referred to as almost (or quasi-) incompressibility whenever the term Ψ vol yields J ≈ 1. More details can be found in [Hol02].

Numerical solvers
In this section we briefly review the nonlinear methods we consider, namely Newton-Krylov and quasi-Newton methods. Details regarding their geometric motivation can be found in [WN99], whereas the convergence properties are detailed in Section 5. Our study is mainly motivated by the use of these methods in an HPC infrastructure, so the use of iterative solvers within the nonlinear method is not only inevitable, but also desirable due to their well-established parallel performance. Inspired by this, we have classified our methods according to the degree of exactness of the linear solver it uses, so that a fully inexact method will simply consider the action of the associated preconditioner instead of solving a linear system. Instead, a quasi-exact 1 method will consider an accurate linear solver given by small tolerances, i.e. 10 −10 and 10 −8 absolute and relative tolerances, respectively. Now we provide more details regarding the methods, and in particular emphasize where the linear solver is used.
Newton-Krylov methods. These methods can be seen as either a minimization procedure or as a root finding algorithm, the latter often referred to as Newton-Raphson method. We present it as a root finding method. Consider Equation (1) in quasi-static, residual form with P = ∂Ψ ∂F for a given Helmholtz potential Ψ, and consider an initial iterate d 0 . Then, given a previous iteration d n−1 , the next one is given as the solution δd n of the linearized problem where ∂ d stands for the Frechét derivative with respect to d, and the update is then given by d n = d n−1 + δd n . Equation (2) gives rise to a linear system of equations often referred to as the tangent problem or simply the linearized problem, which is solved by means of an iterative Krylov space method, accelerated by a preconditioner, see e.g. [CGKT94, CGK + 98]. In practice, the most used one is referred to as Newton-MG, meaning that problem (2) is solved with an iterative method, preconditioned by an algebraic multigrid preconditioner (AMG) [Hac13,CFPS15].
Inexact Newton-Krylov methods. Problem (2) is an approximation of the actual equation R(d) = 0, and thus it might not be true that an accurate solution of the tangent problem will also yield an accurate solution of the original equation. In fact, this usually gives rise to over-solving the problem [EW94] in the first iterations, which motivates the use of inexact solvers for (2), i.e. to use large tolerances for the solution of the tangent problem. In addition, the quality of the linearization improves as the iterates are closer to the solution, which motivates the use of adaptive (relative) tolerances. In particular, we use the Einstat-Walker strategy [EW96], given by The choice of the norm · is arbitrary, so we consider in what follows the 2 norm. This of course makes sense only in the discrete setting, meaning that for each residual vector r we consider the norm r = i r 2 i 1/2 . The resulting scheme guarantees superlinear convergence [DES82], which is of course worse than the quadratic convergence of a standard Newton-Krylov scheme, but gives an overall reduced complexity as it avoids oversolving the linearized problem.
BFGS method. This is a Quasi-Newton method, and although it was initially devised as a minimization procedure [WN99], it can also be adapted to be used as a root-finding algorithms [DS96]. For this, we look at the minimization principle from which (1) in quasi-static form is obtained: A quadratic approximation of this problem yields and such scheme is indeed equivalent to (2) when applied as an iterative procedure with an exact Hessian. This method requires an initial approximation of the Hessian B −1,0 ≈ [∂ 2 d Π(d 0 )] −1 , which is then enriched at each iteration with a rank two perturbation that includes curvature information given by The approximate Hessian is then used to compute the next iteration as in Equation (2): The action of B −1,k+1 is implemented by means of a two-level recursion that allows for a limited memory implementation [Noc80]. This method, as the inexact Newton-Krylov, yields superlinear convergence, and has the additional advantage of requiring only one assembly of the Hessian matrix (or any other initial matrix). The initial approximation of the Hessian is critical for the convergence and performance of this method, and we have indeed observed that simpler approaches do not yield satisfactory results for this method in the context of nonlinear mechanics. Motivated by this, we leverage the preconditioners obtained from the initial matrix ∂ 2 d Π(d 0 ) to obtain better approximations of the Hessian. If we denote by P the preconditioner arising from the initial Hessian, we consider the following approximations: • B −1,0 = action of P obtained from Hessian matrix (BFGS-preonly) • B −1,0 = inexact Krylov solver, large (≈ 10 −2 ) relative tolerance (inexact BFGS) • B −1,0 = quasi-exact Krylov solver, small (≈ 10 −6 ) relative tolerance (quasi-exact BFGS) We have added on the right the names corresponding to Figure 1. We note that the quasi-exact and inexact versions are not covered by the theory, as the action of B −1,0 changes at each iteration. This happens because the number of linear iterations required at each nonlinear iteration varies according to the tolerance. We highlight that the use of a fixed tolerance in the inexact scenario is considered for simplicity, as an equivalent Einstat-Walker type of adaptive tolerance could be considered. This is an interesting alternative and by no means a trivial one to set up, so we leave this for future work.
We show a classification of these methods according to the exactness of the linear solver in Figure 1, where we have fixed three degrees of exactness: preconditioner only (preonly 2 ), inexact and quasi-exact. The two families present different types of performance: the Newton family relies on a low iteration count, where the level of inexactness relaxes the linear solver at each iteration, at the cost of some additional Newton steps. The quasi-Newton family instead yields a larger number of iterations, but each iteration is much cheaper than a Newton iteration, as it does not requires the Jacobian matrix to be reassembled. The level of exactness in this case reduces the iteration count, at the cost of making the iterations more expensive. The extreme cases Newton-preonly and quasi-exact BFGS have shown to be noncompetitive in our experiments, so we do not consider them in what follows. Throughout the remaining parts of this work, whenever no confusion arises, we may refer to Newton-Krylov and inexact Newton-Krylov methods as simply Newton methods. Throughout this manuscript, we will denote the Newton-Krylov, inexact Newton-Krylov, BFGS-preonly (or simply BFGS) and inexact-BFGS methods as NK, iNK, B and iB respectively, as shown in Figure 1. We also show in Appendix A the commands used to run each of the four methods under consideration.
Less inexact More inexact smaller ← Relative tolerance of linear solver → larger Figure 1: Description of exactness level in the linear solver and nonlinear solvers used in the tests, where we denote "preconditioner only" as "preonly". Newton-preonly and quasi-exact BFGS have been discarded from our experiments because they are not competitive: the first requires too many matrix assembles, the second spends too much time at each iteration.

Test problems and solvers details
In this section, we describe the test cases we use to compare the performance of each method. The scope of the heartbeat test is that of providing a more realistic scenario in which the insight acquired in the other two tests will be transfered. Details on the implementation will be given on each problem.
Cook test. This is a static almost incompressible mechanics benchmark [PM16]. We use homogeneous Dirichlet conditions on {x = 0}, a vertical traction given by t = (0, τ, 0) with τ = 10 6 in deformed configuration on {x = L}, a homogeneous distributed load f = 0, and homogeneous Neumann conditions on the remaining parts of the boundary. The constitutive modeling is given by a almost incompressible Neo-Hookean material: where C 1 = 1 2 µ, µ = 8.194 · 10 7 , k = λ + 2 3 µ, λ = 2µν/(1 − 2ν), and ν = 0.3. The Krylov solver was configured with a GMRES linear solver without restart and with a modified Gram-Schmidt procedure which is more robust [GVL96], preconditioned with HYPRE-BoomerAMG [FY02], using its default configuration and setting the matrix block size to 3 in PETSc to improve the efficiency of the preconditioner. The absolute tolerance was set to 10 −8 for the exact solver, and to 0 for the inexact ones, whereas the relative tolerance was set to 10 −6 for the exact solver and to 10 −2 for the inexact ones. We show the solution in Figure 2 (a).
Twist test. This is a static incompressible mechanics problem [BGO15], where we use the inf-sup stable stable finite elements P 2 − P 0 for the approximation of the displacement and pressure. Boundary conditions are given by a homogeneous on {z = 0}, a π/6 rotation in the x and y components on {z = L z }, homogeneous data f = t = 0, and homogeneous Neumann conditions on the remaining parts of the boundary. The constitutive modelling is given by the following polyconvex potential: where α p = β p = α s = β s = 9000 and cof F := det(F)F −T . The logarithmic term with det(F) is required because of the weak imposition of the incompressibility constraint [ADVL + 13]. The preconditioning of a saddle point problem is more challenging, so to obtain a more robust performance we used a lower Schur complement preconditioner based on a (d, p) field split [BKM + 12], we provide details on the Schur complement preconditioner in Appendix B. The displacement block uses an AMG preconditioner, the Schur complement block instead uses the SIMPLE preconditioner from fluid mechanics [EHS + 06], meaning that the inverse of the displacement block is approximated by the inverse of its diagonal, and for the resulting block we used a block Jacobi preconditioner, which we observed to suffice. Two additional comments are in place: (i) to calibrate the effectiveness of each solver in the Schur complement, we started from using a direct solver in each block to guarantee convergence in at most two iterations [Man90], and started relaxing the blocks from there. The difficulty of the problem is predominantly the displacement block, which can be appreciated by the simple preconditioner used for the Schur complement block; (ii) the pressure block presents much less degrees of freedom (DoFs), so we used a telescopic approach that creates an MPI subcommunicator with a fixed 25% of the original processes to avoid excessive communication. We show the solution in Figure 2 Heartbeat test. This problem represents a real-case study and is indeed our case of interest. It models the contraction of a human left ventricle in an idealized geometry given by a prolate ellipsoid, with a pointwise set of coordinates representing the muscle fiber orientation (f 0 , s 0 , t 0 ) obtained through rule based methods [BBPT12]. We consider a Guccione hyperelastic potential [GMW91], given by which are the components of E in the fiber-induced frame of reference {f , s, n}. See [GMW91] for reference values of the related parameters. The Piola stress tensor is enriched with a fiber-wise force known as active stress that models the contraction of the muscle cells (cardiomyocites): where γ is a given function known as the activation function. For simplicity we consider an analytical activation given by where the peak activation constant is C PA = 10 4 and the period is T = 0.8. More details on the generation of the fibers and the activation function can be found in [BBPT12] and [RSA + 20] respectively. Boundary conditions are given by homogeneous Neumann on the endocardium ∂Ω endo and the base ∂Ω base , whereas the epicardium ∂Ω epi considers a Robin condition that models the friction of the epicardium with the pericardium [ULM02] and is given by More information of the physical meaning of these boundaries and the motivation for these choices can be found in [QLRRB17]. The Krylov solver in this case is identical to the one used in the Cook test. We show the solution at t = 0.2 in Figure 2 (c).
(c) Heartbeat test. Figure 2: Solutions of the benchmark tests used. Color is scaled according to displacement norm |d|, and the reference configuration is depicted with its mesh. In (c), the mesh is cut roughly at half the height of the geometry.

Theoretical results and solvers expected behavior
The purpose of this section is to show that the hypotheses from the convergence analysis, despite not being directly applicable to this problem, provide fundamental insight to guide the choice of the nonlinear solver.
To state the results we consider an arbitrary scalar function f : R m → R to be minimized and an initial point x 0 . The results and their proofs can be found in [WN99]. First we state the convergence result of Newton's method, where we remark that the inexact variant requires no additional hypotheses.
Theorem 1 (Newton methods) Assume that there exists a point x * such that ∇f (x * ) = 0, ∇ 2 f is Lipschitz continuous in a neighborhood of x * and ∇ 2 f (x * ) is positive definite. Then, x * is a local minimizer and, if the initial guess x 0 is sufficiently close to x * , the iterates obtained by Newton's method converge quadratically to x * . Now we state the general convergence theorem for the BFGS method.
Theorem 2 (Quasi-Newton methods) Assume that f is twice continuously differentiable and that ∇ 2 f is Lipschitz continuous at the minimizer x * . Assume also that the level is convex, that ∇ 2 f is positive definite on L and that the initial matrix B 0 is positive definite. Then, the sequence of iterates obtained by the BFGS method converges superlinearly to x * .
We make the following observations, which are fundamental to understand the performance of these methods.
• Both methods require the initial point to be sufficiently close. This can be seen in the BFGS method through the convexity of L.
• The BFGS method is better suited for convex problems, given that the Hessian is required to be definite positive not only at the minimizer as Newton, but also in the entire level set L.
• Inexact Newton-Krylov convergence holds under the same hypotheses of Newton's method, plus some hypotheses on the relative tolerances η k that guarantee the superlinear convergence. We consider tolerances that satisfy such hypotheses [EW94], and note that they can sometimes yield global convergence, unlike Newton's method. This makes inexact Newton-Krylov methods attractive as they can be potentially more robust than a classic Newton method.
• The initial matrix approximation B 0 is variable in the inexact BFGS method, meaning that its performance might not necessarily be better, nor more robust, than the case in which only the action of the preconditioner is considered. This is due to this case not being covered by the theory, so that in principle there could be additional hypotheses on the action of a variable B 0 .
• Notes for different types of problems: -Cook test: In this case, the problem being solved is polyconvex, but if a sufficiently close initial guess is considered, we can possibly fall in an attraction basin. This means that both Newton and BFGS methods should convergence, but both BFGS might possibly be less robust due to the convexity requirement of L.
-Twist test: This is a saddle point problem, so BFGS is not guaranteed to converge even if a sufficiently good initial approximation is considered. This can be seen from the positive definiteness hypothesis of the Hessian on L.
-Heartbeat test: The inertia term in elastodynamics results in a convex contribution to the variational principle associated to the problem, so we expect this case to be easier to solve in spite of the more complex nonlinearities involved. Instead, the nonlinearities should impact the number of nonlinear iterations.

Numerical results
In this section, we compare the performance of the solvers introduced in Sec. 3, for all tests under consideration. The scope of these tests is to assess the sensitivity of the methods with respect to the problem size, some of the problem parameters, and the number of computational cores. For the heartbeat problem, we consider four different meshes with the degrees of freedom shown in Table 1. We note that meshes 3 and 4 are refinements of meshes 1 and 2 respectively. Note as well that since BFGS-preonly does not solve any linear system, we do not report the linear iterations for this method.

Sensitivity with respect to problem size
To compare the performance of all methods, we report the total nonlinear iterations counts, the average Krylov iterations per nonlinear step and the solution time as we increase the number of degrees of freedom.
Cook test. We present the nonlinear iterations, linear iterations and CPU time in Table 2, varying the dimension of the problem from about 10 thousand to 1.5 million degrees of freedom. We first observe that notably all methods exhibit a robust behavior with respect to the degrees of freedom, except for the BFGS-preonly with second order finite elements. Indeed, the total and average linear iterations reported are roughly constant, but in general the inexact-BFGS yields lower average linear iterations. All methods exhibit an approximately linear increase of CPU time with respect to the degrees of    freedom, with the exception of BFGS-preonly with second order elements. In general, BFGS methods outperform Newton methods in terms of the execution time, but attention must be given to the second order scenario. The ratio between the CPU times of the best method (BFGS-preonly) and Newton's method in the case with the largest amount of DoFs is 0.63 for P 1 and 0.67 for P 2 elements.
Twist test. We report the nonlinear iterations, average linear iterations per nonlinear step, and CPU times in Table 3. We observe that BFGS methods do not converge for this problem, so we can not report results for them. This holds not only for both methods presented (BFGS, inexact-BFGS), as we also tried the BFGS-exact variant described in Section 3, which was ineffective as well. Both Newton-Krylov and inexact Newton-Krylov methods present a very mild increase in the nonlinear iterations, but the inexact variant presents a plateau of average linear iterations, whereas the Newton-Krylov method presents a monotonic increase in the average linear iterations. This results in the inexact method being a faster method, and the difference in CPU times increases with the problem size. The ratio between the CPU times of the inexact Newton-Krylov method and Newton-Krylov method in the case with the largest amount of DoFs is 0.79.
Heartbeat test. We report the nonlinear iterations, average linear iterations per nonlinear iteration, and CPU time for first and second order finite elements averaged during the first timestep in Table 4 and their time evolution in Figure 3. We note that all methods are robust with respect to the degrees of freedom, except the inexact-BFGS, as it does not converge for fine meshes when using second order elements. Both Newton methods present similar nonlinear iteration counts, and all methods have a similar average linear iteration count. This results in the inexact-BFGS having a lower average linear iteration count overall as it requires more nonlinear iterations. This has an impact in the overall solution time, where BFGS methods perform much better than Newton methods. The ratio between the CPU times of the best method (inexact-BFGS in P 1 , BFGS-preonly in P 2 ) and the Newton-Krylov method is 0.63 for P 1 and 0.57 for P 2 elements.

Robustness with respect to the data
To study the robustness, we varied in each test case the parameter which resulted in a larger deformation. In the Cook test this refers to the vertical load τ , in the Twist case this refers to the angle of rotation in the boundary conditions, and on the Heartbeat case this refers to the peak activation constant C AP .
Cook test. We vary the magnitude of the distributed traction term τ on the right from 1.5 · 10 6 Pa to 2 · 10 6 Pa, where values over the upper bound considered make all methods diverge. The nonlinear iterations, total linear iterations and average linear iterations per nonlinear step are shown in Figure  4. We have computed the solution for four values of τ , and plotted the result only when the solver converged. Most notably, only the inexact Newton-Krylov method is able to solve the highest load value, confirming that it can show more robust behavior than the Newton-Krylov method. Moreover, the inexact-BFGS does not converge for some values in case of second order elements, making it the least robust method in this scenario. In general, the difference in robustness is much clearer with second order elements, as also BFGS-preonly and Newton-Krylov converge only for two of the four values considered for τ . When looking at the linear iteration counts, we see that the methods that converge tend to maintain a roughly constant number of average linear iterations, so that the increase in total linear iterations is given by the additional nonlinear iterations incurred. The only exception is the Newton-Krylov method for second order elements. Interestingly, inexact-BFGS is more robust than BFGS-preonly for first order elements, but instead BFGS-preonly is slightly more robust than inexact-BFGS for second order elements. There is no advantage in considering a Newton-Krylov method, but instead its inexact variant is much more robust with respect to τ .
Twist test. We vary the angle in the boundary conditions from π/6 to π/2, and show the results in Figure  5. We note that both Newton methods converge for all values under consideration, and diverge when higher ones are considered, so there is no significant difference in terms of robustness between them. Still, the Newton-Krylov method presents and increase in the total linear iterations, whereas the inexact variant shows a more robust iteration count.
Heartbeat test. We vary the peak activation C P A from 10 4 to 10 5 , and plot only the results when the methods converged. In this case, we note that Newton methods provide an improved robustness for both first and second order finite elements, and as in the other tests, the inexact variant is more robust. This effect is inverted in the BFGS methods, which are less robust and instead the BFGSpreonly shows an improved robustness with respect to the activation force. We again observe that the total linear iterations of the Newton-Krylov method increase with the problem difficulty, whereas the inexact Newton-Krylov method maintains a roughly constant number of linear iterations despite the increase of the problem difficulty.

Scalability
We test the strong scalability with respect to the finest mesh used in the performance tests. In the heartbeat test, we consider only the first 10 timesteps.
Cook test. We report the nonlinear and average linear iterations and the CPU times in Table 5. The CPU times are also displayed in Figure 7. We note that all methods are robust with respect to the number of cores, as could have been expected in virtue of the robustness of the preconditioner. Still, it is very interesting to see that there is a point in which BFGS-preonly converges with a much lower iteration number in the P 1 case, namely when using 8 processors. This is the only case in which the CPU time is comparable to that of the inexact-BFGS method. For any other number of cores, the inexact-BFGS outperforms all other methods for both P 1 and P 2 elements. As already observed in the performance test, the BFGS-preonly method is less robust than the other methods.
Twist test. We report the nonlinear and average linear iteration counts, together with the CPU times in Table 6. The CPU times are also displayed in Figure 8. We note that this test in particular is in complete agreement with the previous ones, meaning that both Newton-Krylov methods perform similarly and exhibit an overall robust behavior. Both methods present adequate strong scaling, with the inexact method showing better solution times.
Heartbeat test. For this test, we report the nonlinear and average linear iterations, and CPU times averaged for the first 10 timesteps in Table 7. The total CPU times as function of the number of cores are displayed in Figure 9. We can appreciate that all methods are robust with respect to the number of cores, as these yield no significant variations. From the scalability curves shown in Figure 9, we note that all methods show adequate scaling, with the BFGS methods presenting an overall better performance. More specifically, inexact-BFGS is faster for P 1 , and BFGS-preonly is faster for P 2 .

Conclusions
This work presents a detailed numerical study of some of the main numerical solvers used in large scale simulations of nonlinear mechanics, where all of the fundamental components of the methods are included: the nonlinear solver, the linear solver and the preconditioner. The main finding from this study is that superlinear solvers such as inexact Newton-Krylov and BFGS methods present an overall robustness similar to the gold standard represented by the Newton-Krylov method in nonlinear elasticity, but with improved solution times.
In difficult problems such as static and incompressible mechanics, the choice of the solver can be defined a-priori according to the convergence hypotheses required for each method, where in fact we have seen the BFGS methods fail for incompressible mechanics, as they usually require some form of convexity to guarantee convergence. Indeed, BFGS methods showed the best solution times in static mechanics, while inexact Newton-Krylov showed the best performance for the incompressible tests. The time-dependent heartbeat test showed us that the convexity contributed by the inertia term can further improve the performance of the BFGS method, which showed much more robust iteration counts in this scenario than those shown in the static mechanics case. This comes as a surprise, given that the constitutive model used in the heartbeat test is much more complex than the models of the other tests. These nonlinear solvers present no drawbacks in terms of both performance and strong scalability, suggesting that in the context of nonlinear elasticity inexact Newton-Krylov and quasi-Newton methods should be preferred to standard Newton-Krylov methods.
We conclude by mentioning two possible extensions of this work. The first one is the inclusion in this study of Jacobian lagging techniques [BB13], which consist in avoiding the reassembly of the Jacobian for some Newton iterations. These lagging techniques could potentially further improve the performance of the Newton methods considered. Instead, inspired by our inexactness taxonomy, we note that a Jacobian reassembly could be considered for the BFGS methods, where the Jacobian gets sometimes reassembled to improve the quality of the initial Jacobian. Both strategies are non trivial to implement, and require a detailed problem-specific study to obtain the best performance. The second interesting extension is the inclusion in our nonlinear mechanics study of a polyconvex potential [Bal76], even if solvers for this specific type of potential are not yet available, to the best of our knowledge. One interesting work in this direction is [BGO15], where a Hu-Washizu [Was68] formulation is used to exploit the convexity of the auxiliary function arising from the polyconvex potential. Still, this formulation results in a problem with many auxiliary variables, and its efficient numerical approximation remains a research topic.

A PETSc options
The commands used to invoke each of the methods described in Section 3 is detailed in what follows.
• Inexact Newton-Krylov (iNK) -snes_type newtonls -ksp_type gmres -pc_type hypre -ksp_atol 1e -14 -ksp_rtol 1e -2 # Should not be used , set for safety -snes_ksp_ew -sn e s _ k s p _ e w _ r t o l 0 1e -1 -s n e s _ k s p _ e w _ r t o l m a x 0.1 -s n e s _ l i n e s e a r c h _ t y p e basic Listing 2: PETSc commands to use inexact Newton-Krylov.
• BFGS-preonly (B) -snes_type qn -ksp_type preonly -pc_type hypre -snes_qn_type lbfgs -snes_qn_m 50 -s n e s _ q n _ s c a l e _ t y p e jacobian -s n e s _ l i n e s e a r c h _ t y p e basic Listing 4: PETSc commands to use BFGS.

B Schur complement preconditioners
A Schur complement preconditioner is one that arises from a block LU factorization according to two (arbitrary) index sets of a matrix. If we consider a block matrix M given by the general structure with A invertible, a Gaussian elimination procedure yields We note that if C is invertible, the same procedure can be applied with respect to it. Schur complement based preconditioners enjoy excellent theoretical properties, as the preconditioned system possesses at most 3 distinct eigenvalues [MGW00], implying that it converges in at most 3 iterations of a Krylov subspace method. This is true whenever the Schur complement S = C−B 2 A −1 B 1 is evaluated exactly, which is usually computationally intractable. One simple approximation, which is the one we use for the preconditioning the incompressible mechanics problem, is to consider the approximation A −1 ≈ diag −1 (A). Another important point is that in general using all three blocks arising from the factorization in (5) is not necessary, and instead it is sufficient to consider a lower factorization (first two blocks) or an upper factorization (last two blocks). Whenever C = 0, it is also possible to use a diagonal factorization (middle block only).