Truncated nonsmooth Newton multigrid for phase-field brittle-fracture problems, with analysis

We propose the truncated nonsmooth Newton multigrid method (TNNMG) as a solver for the spatial problems of the small-strain brittle-fracture phase-field equations. TNNMG is a nonsmooth multigrid method that can solve biconvex, block-separably nonsmooth minimization problems with linear time complexity. It exploits the variational structure inherent in the problem, and handles the pointwise irreversibility constraint on the damage variable directly, without regularization or the introduction of a local history field. In the paper we introduce the method and show how it can be applied to several established models of phase-field brittle fracture. We then prove convergence of the solver to a solution of the nonsmooth Euler–Lagrange equations of the spatial problem for any load and initial iterate. On the way, we show several crucial convexity and regularity properties of the models considered here. Numerical comparisons to an operator-splitting algorithm show a considerable speed increase, without loss of robustness.


Introduction
The equations of phase-field models of brittle fracture present a number of challenges to the designers of numerical solution algorithms [2].Even in the smallstrain case the equations are nonlinear, due to the multiplicative coupling of the mechanical stresses to the degradation function.At the same time the non-healing condition introduces an inequality constraint.Finally, eigenvalue-based splittings of the energy density as in [36] make the equations nondifferentiable.
In this paper we focus on the spatial problems of small-strain brittle-fracture phase-field models obtained by a suitable time discretization.The standard approach to solving these spatial problems is based on operator splitting.Algorithms based on this approach, also known as staggered schemes, alternate between solving a displacement problem with fixed damage and a damage problem with fixed displacement.Both subproblems are elliptic and well-understood, and such methods are therefore straightforward to implement.The method can be interpreted as a nonlinear Gauß-Seidel method [14], which provides a natural framework for convergence proofs.Applications of the operator splitting scheme and its extensions appear, e.g., in [7,9].With particular semi-implicit time discretizations it is also possible to solve the spatial problems by solving only one damage problem and one The authors gratefully acknowledge the financial support by the German Federal Ministry of Education and Research through the ParaPhase project within the framework "IKT 2020 -Forschung für Innovationen" (project numbers 01-H15005C, 01-15005D, 01-15005E).
displacement problem [35].This is very fast, but works only if the load steps are small enough.
In contrast, other works propose monolithic solution schemes based on Newton's method [14,15,49,51,52].For the unmodified Newton method only local convergence can be shown, and failure to converge for large load steps is readily observed in practice [54].Therefore, various authors have proposed extensions or modifications of the Newton idea to stabilize the method.In [15] a line search strategy is applied to enlarge the domain of convergence of Newton's method.Wu et al. [54] and Kristensen and Martínez-Pañeda [27] propose to use the BFGS quasi-Newton algorithm, claiming that it is more stable than Newton's method and more efficient than operator splitting.Heister et al. [22] proposed a modified Newton scheme which was later improved by Wick [52] with an adaptive transition from Newton's method to the modified Newton scheme.Recent results in [25] suggest that nonlinear preconditioning can speed up the solution process tremendously.Finally, the authors of [34,45] suggest an arc-length method based on the fracture surface, and an adaptive time stepping scheme to enhance the robustness.In summary, while monolithic Newton-type methods are reported to be faster than operator-splitting algorithms, the latter ones are more robust [2,54].
Various approaches are used in the literature to deal with the damage irreversibility.A natural approach is to regularize the constraint, as investigated in [25,36,50].This leads to an additional parameter, and to ill-conditioned tangent matrices [16].Interior-point solvers implement an automatic control of the new parameter; they are investigated in [49].An alternative approach considers the thermodynamic driving force of the fracture phase-field as a global unknown yielding a three-field formulation which results in a saddle-point principle [36].A third formulation considers the Karush-Kuhn-Tucker conditions and shifts the thermodynamic driving force of the fracture phase-field into a local history field representing the maximum over time of the elastic energy [35].This approach, frequently known as the Hfield technique, therefore trades the inequality constraint for the nondifferentiable maximum function.Alternatively, the time discrete problems can be reformulated as semismooth systems by means of so-called complementarity functions.This strategy can be combined with monolithic semismooth Newton techniques [32] or nested Newton and active set methods [22].Unfortunately, these approaches spoil the variational structure of the spatial problems.Augmented Lagrangian solvers as in [50] introduce extra variables.Closest in spirit to the present manuscript is the use of bound-constrained optimization solvers, used, e.g., in [3,9,14,53].None of these approaches are fully satisfactory.
The effect of the nondifferentiable terms caused by anisotropic splittings of the mechanical energy density as in [3,36,46] is rarely discussed in the literature.Hybrid formulations like the one proposed in [2] try to overcome the additional computational difficulties of these splittings by further changes to the model, again at the cost of sacrificing the variational structure.
All these approaches are slow in the sense that they have to solve global partial differential equations at each Newton or operator-splitting iteration.This is expensive, even if efficient multigrid methods are used for the linear subproblems (as, e.g., in [23]).When the methods use direct sparse solvers for the linear tangent problems, memory consumption can become problematic, too.At the same time, the problem of small-strain phase-field brittle fracture has a lot of elegant variational structure; in particular, it fits directly into the rate-independent framework of Mielke and Roubíček [37].As a consequence, implicit time discretization leads to a sequence of coercive minimization problems for the displacement and damage fields together.These problems are not convex, but they are biconvex, i.e., convex (even strongly convex) in each variable separately.Pointwise inequality restrictions ḋ ≥ 0 to handle the irreversibility of the fracture process as proposed in [36] reduce the smoothness of the objective functional, but do not influence its convexity or coercivity properties.The same holds for anisotropic energy splittings based on linear quantities or the eigenvalues of the mechanical strain.
Recent years have shown that nonsmooth multigrid methods are able to solve variational nonsmooth problems from mechanics efficiently without the need for solving global linear systems of equations [18,20,21,24,26,44].This can make them vastly more efficient than operator-splitting or Newton-based methods.As there are no sparse matrix factorizations, memory consumption remains linear in the number of unknowns.In addition, these multigrid methods can be shown to converge globally (i.e., from any initial iterate and for any load step) to a stationary point of the objective functional.The proof exploits the above-mentioned variational structure, together with certain separability properties.As one such method, the Truncated Nonsmooth Newton Multigrid method (TNNMG) can treat the pointwise constraints of the increment problems directly, i.e., without artificial regularization or tricks like the H-field technique [35].The idea is that TNNMG only needs to handle these constraints in a series of low-dimensional subproblems, each of which is easy to solve by itself.As a consequence, solving the problems with constraints is not appreciably slower than solving the corresponding unconstrained problem.
In this paper we show how the TNNMG method can be used to solve small-strain brittle-fracture problems.This involves in particular verifying that the increment functionals have the required convexity and smoothness properties.We do this for a range of different degradation functions and local crack surface densities (including the standard Ambrosio-Tortorelli functionals of type 1 and 2), closely related to the family of models considered in [10].We cover elastic energies with various types of anisotropic splittings, including the splitting based on strain eigenvalues of Miehe et al. [36].For the proofs we use results from the convex analysis of spectral functions [4,41].Extension to the slightly more general damage models of [33] is straightforward, provided the stored elastic energy has the required convexity and smoothness properties.In contrast to the multilevel trust region method proposed in [24], the TNNMG method presented here relies on nonsmooth Newton techniques leading to linear subproblems and thus gives more flexibility in the selection of coarse grid solvers.
The paper is organized as follows: Section 2 discusses a framework of small-strain phase-field models for brittle fracture, and shows the range of applicability of the TNNMG solver.Section 3 introduces the natural fully implicit time discretization, and proves existence of solutions for the spatial problems.In both sections we pay particular attention to the mathematical properties of the energy functionals.In Section 4, finally, we introduce the TNNMG method.We explain its construction, discuss various algorithmic options, and prove that it converges globally to stationary points of the increment energy functional.The numerical efficiency is then demonstrated in Section 5. Our reference for comparison, briefly revisited in Section 5.1, is the operator-splitting iteration proposed in [10], which uses a projected Newton method [6] for the constrained damage problems.We compare the solvers for two-and three-dimensional example problems with different forms of the local crack surface density, and with and without spectral splittings.We observe a noticeable performance increase, without loss of robustness.

Phase-field models of brittle fracture
This section presents a range of phase-field models for brittle fracture, and discusses its smoothness and convexity properties.
Consider a deformable m-dimensional object represented by a domain Ω ∈ R m .The deformation of such an object is characterized by a displacement field u : Ω → R m .The object is supposed to exhibit small-strain deformations and elastic material behavior only, and we therefore introduce the linearized strain tensor ε(u) := 1 2 (∇u + ∇u T ).Following [36], we model the fracturing by a scalar damage field d : Ω → [0, 1], where d = 0 signifies intact material, and d = 1 a fully broken one.Dirichlet boundary conditions can be posed both for the displacement and for the damage field.For this we select two not necessarily equal subsets Γ D,u , Γ D,d ⊂ ∂Ω of the domain boundary, and require where u 0 and d 0 are two given functions.Displacement and damage field evolve together, governed by a system of coupled nonsmooth partial differential equations.Disregarding inertia effects, we obtain a rate-independent system in the sense of Mielke and Roubíček [37].Such a system can be written using the Biot equation In this equation, E is a potential energy, which we assume to be of the form The term ψ is a degraded elastic energy density, and will be discussed in detail in Section 2.1.The term γ models the local crack surface density, and will be discussed in Section 2.2.The number g c is Griffith's critical energy release rate, a material parameter.P ext represents time-dependent volume and surface forces, which drive the evolution.We assume that P ext is linear and H 1 (Ω)-continuous in u, and differentiable in t with bounded time derivative.The last term of (2) implements the restriction that the damage field can only assume values between 0 and 1.For a set K ⊂ R we define the indicator functional For a closed, convex, nonempty set K, the functional I K is convex, lower semicontinuous, and proper.Adding the constraint d ∈ [0, 1] explicitly is not always necessary, as some fracture models lead to evolutions that satisfy the constraints implicitly.However, as pointwise bounds come with practically no cost when using the TNNMG solver, we do include them to extend our range of models.
To make the potential energy E well defined, we will in general consider it on the first-order Sobolev space H 1 (Ω, R m ) × H 1 (Ω, R).Incorporating the boundary conditions leads to the affine subspace The second term of the Biot equation ( 1) is It implements the pointwise non-healing condition ḋ ≥ 0 as proposed by [36].Note that R(d, •) : H 1 (Ω) → [0, ∞] is convex and lower semicontinuous, and R(d, 0) = 0.The fact that R is positively 1-homogeneous in ḋ implies the rate-independence of the system.Since the particular functional (3) only depends on ḋ but not on d, we will also write R( ḋ) := R(d, ḋ) and ∂R( ḋ) := ∂ ḋR(d, ḋ).
Remark 2.1.Using the definition of the subdifferential ∂R it is straightforward to see that the Biot equation ( 1) is equivalent to the variational inequality Likewise, it is equivalent to the coupled system ∀e ∈ H 1 0 , where we have denoted by H 1 0 ×H 1 0 := H 1 u0 ×H 1 d0 −(u 0 , d 0 ) the homogeneous space corresponding to H 1 u0 ×H 1 d0 .Since these two notions of solutions are based on firstorder derivatives of E and R, they may lead to local minimizers or saddle points of the functional E. To overcome this, there are alternative so-called energetic formulations of the problem.In general energetic solutions are solutions of (4), while the converse is only true for convex E. We will not discuss such solutions here, but refer to [37,47], where they are discussed for damage-related and more general rate-independent processes.Remark 2.2.In the engineering literature, the same problem is sometimes formulated as Requiring this to contain 0 is the Biot equation (1).
We now discuss the potential energy and the dissipation potential.
2.1.Degraded elastic energy density.We consider models that behave linearly elastic and isotropic if the material is in an undamaged state.That is, for the undamaged stored energy density we use the St. Venant-Kirchhoff material law, whose energy density is given by with Lamé parameters µ > 0 and λ > − 2 3 µ.With this choice of parameters the quadratic functional ψ 0 is strongly convex on S m , the set of real symmetric m × m matrices.
The undamaged energy density is splitting as into a part ψ + 0 that produces damage and another part ψ − 0 that does not.The damage-producing part is then scaled by a so-called degradation function and the energy density ψ : The residual stiffness k > 0 guarantees a well-posed problem in case of fracture.
Remark 2.4.As an alternative to this assumption, one could consider degradation functions with g(0) = 1 − k such that the original energy density ψ 0 is exactly recovered in the undamaged case (see, e.g., [51]).However, in the typical regime of k 1 both will essentially yield the same result.Thus we will follow the more common approach g(0) = 1 here, although the proposed methods could equally be applied to the alternative one.
Note that several authors require g (1) = 0 in order to ensure that the evolution does not lead to values of d larger than 1.We do not need this assumption here, because the pointwise constraint d ≤ 1 is enforced explicitly by the energy term (2).
Note that the functions g a and g d are strictly convex, but g b and g c are not even convex.For the rest of the paper we will restrict our considerations to convex twice continuously differentiable degradation functions g.Various splittings of ψ 0 have been proposed in the literature.We cover four common strain-based splittings taking the form (5). 1 All those splittings have the property that ψ 0 = ψ + 0 + ψ − 0 , and we will show that all have the following essential properties:  [42]).Notice that the strong convexity (P 4 ) implies strong monotonicity of ∇ψ, i.e., Here and in the following we denote by •, • : V * × V → R the duality pairing of a vector space V. 1 The stress-based splitting of Steinke and Kaliske [46] is left for future work. 2 While we only require L ≥ 0 in (P 3 ), the strong convexity (P 4 ) in fact implies L > 0.
For the splittings considered in the following we will only prove (P 1 ) and (P 2 ) directly, and show that the simplified assumptions of the following lemma hold true.This then implies (P 3 ), (P 4 ), and (P 5 ).
To show strong convexity, we first note that ψ 0 = ψ + 0 + ψ − 0 is strongly convex on S m with a modulus η > 0 independent of d.Now consider the function Since this is a weighted sum of two convex functions ψ + 0 and ψ − 0 with non-negative weights g(d) ≥ 0 and 1, it is itself convex.Thus, as a sum of this convex function and the strongly convex functions Cψ 0 , the function ψ(•, d) is itself strongly convex and inherits the convexity modulus kη of kψ 0 .Finally, we note that with the same η we have Despite those strong properties of ψ(•, d) we note that ψ(ε, d) is not convex in d and ε together for any of the splittings considered below.
2.1.2.Volumetric decompositions.The isotropic model is unphysical, because it produces fracturing for all kinds of strain.In [29], Lancioni and Royer-Carfagni obtained better results by letting only the deviatoric strain contribute to the degradation.They introduced the split With these definitions, the energies are Lemma 2.7.The energy density ψ defined in (5) with the isotropic volumetric splitting (7) has the properties (P 1 )-(P 5 ).Furthermore ψ(•, d) has the stronger property that it is in C ∞ and quadratic for all d ∈ [0, 1].
The decomposition of Lancioni and Royer-Carfagni is still isotropic.Amor et al. [3] proposed to only degrade the expansive part of the volumetric strain.Using the ramp functions x + := max{0, x}, x − := min{0, x} that provide the decompositions x = x + + x − and x 2 = x 2 + + x 2 − , they proposed the energy split where only the tensile volumetric strain contributes to damage.
For g(d) + k = 1 the functional ψ(•, d) is quadratic and thus C 2 .In the case g(d) = 1, if ψ(•, d) would be C 2 , then the function t → ψ(tI, d) would also be C 2 .However, this function takes the form and is thus piecewise quadratic but not C 2 in t = 0.
2.1.3.Spectral decomposition.A more elaborate nonlinear splitting separating the tensile and compressive parts of the elastic energy was introduced in [36].To define this splitting it is convenient to introduce the ordered eigenvalue function Eig : S m → R m on the space S m of symmetric m × m matrices, mapping any symmetric matrix M to the vector Eig(M ) ∈ R m containing its eigenvalues in ascending order.Using the ramp functions the tensile and compressive energies ψ + 0 and ψ − 0 are then defined as Note that this indeed defines a splitting ψ 0 = ψ + 0 + ψ − 0 .For this splitting we will make the additional assumption that λ ≥ 0.
To quantify the properties of ψ(ε, d) with respect to the strain tensor ε we use the theory of spectral functions [30].To this end we note that we can write ψ ± 0 as The functions ψ ± 0 are symmetric in the sense that ψ ± 0 (λ) does not depend on the order of the entries of λ ∈ R m .Having this form we can infer properties of the functions ψ ± 0 = ψ ± 0 • Eig from properties of the symmetric functions ψ ± 0 .Lemma 2.9.Let λ ≥ 0. Then the energy density ψ defined in (5) with the spectral splitting (9) has the properties (P 1 )- Proof.We will first show (P 1 )-(P 5 ).An essential ingredient is that the squared ramp functions • 2 ± are non-negative, piecewise quadratic, and convex.(P 1 ) The squared ramp functions • 2 ± and thus Hence the same applies to ψ(•, d).
(  [4] provides global Lipschitz continuity of the gradients ∇ψ ± 0 of the spectral functions ψ ± 0 in the more general context of Euclidean Jordan algebras (which includes the special case of symmetric matrices).In fact, the Lipschitz constant of ∇ ψ ± 0 equals the one for ∇ψ ± 0 if S m is equipped with the Frobenius norm.Using Lemma 2.5 this implies uniform Lipschitz continuity of ψ(•, d).
(P 4 ),(P 5 ) Since the functions ψ ± 0 are weighted sums of convex, non-negative squared ramp functions with nonnegative weights, they are convex and non-negative themselves.Convexity of the functions ψ ± 0 then follows from [4, Theorem 41] while non-negativity of those functions is trivial.Now Lemma 2.5 provides (P 4 ) and (P 5 ).
To characterize second order differentiability of ψ(•, d) we first consider tE, d) for the fixed matrix E with E ij = δ 1i δ 1j would also be C 2 .However, this function takes the form and is thus piecewise quadratic but not C 2 in t = 0.
Remark 2.10.One can show that S m decomposes into finitely many disjoint subsets A i such that ψ(•, d) is twice continuously differentiable in the interior of each of these sets.A matrix ε ∈ S m is in the intersection of several A i if it either has an eigenvalue Eig(ε) i = 0 or if tr ε = 0.While ∇ψ(•, d) is not differentiable at those points, there are still generalized second-order derivatives.For example, the generalized Jacobian in the sense of Clarke contains the derivatives of ∇ψ(•, d) with respect to all the adjacent sets A i .Semismoothness essentially means that such generalized derivatives provide an approximation that can be exploited in a generalized Newton method.
To see this we consider for m = 2 the line segment Then, along this line segment, ψ(•, 1) is quadratic and takes the form which is strictly concave for λ < 0 and sufficiently small k 1.

2.2.
Crack surface density.The crack surface density function per unit volume of the solid is typically of the form [40] γ(d, ∇d) with parameters c w and l, and a parameter function w : [0, 1] → [0, 1].Motivated by the seminal work of Modica and Mortola [38,39], the associated crack surface functional is called a Modica-Mortola functional.The internal length scale parameter l controls the size of the diffusive zone between a completely intact and a completely damaged material.For l → 0 the regularized crack surface yields a sharp crack topology in the sense of Γ-convergence [38,39].For a given function w, the normalization constant c w must be chosen such that the integral of γ(d, ∇d) over the fractured domain converges to the surface measure of the crack set as l → 0. The function w(d) models the local fracture energy.Two types of local crack density functions appear in the literature.Double-well potentials (as briefly reviewed in [2]) provide an energy barrier between broken and unbroken state, but will be disregarded here.Instead, we only consider single-well potentials w, which grow monotonically from the intact state w(0) = 0 to the damaged state w(1) = 1.For such potentials the normalization constant is given by While this follows from the Γ-convergence proof (see, e.g., [1]), the proper constant can also be computed by elementary means: Consider an open set Γ ⊂ R m−1 of measure |Γ| embedded in a domain Ω = R × Γ.We identify Γ with the set {x ∈ Ω : x 1 = 0}, and interpret it as a crack in Ω.To approximate Γ by a damage Using (11) we find that the crack surface energy of the minimizer d is Thus c w has to be selected as in (10) to ensure that the crack surface energy scales like the limit crack surface measure |Γ|.Note that the function d is the w-dependent crack profile, rescaled by the length parameter l.In order to relate the length parameter to the actual crack width, we use the cone construction from [35] and define the crack width to be the average width of a tangential finite cone fitted to the crack profile and the zero function (Figure 1).From the initial condition d(0) = 1, the crack profile equation (11), and the normalization w(1) = 1 it follows that d (0) = −1.Thus-for the normalized solution d-the cone has width 2 at its base and average width 1. Hence the average crack cone width of the rescaled solution d is l.
We will focus on the two widely used potentials and They are referred to in the literature as Ambrosio-Tortorelli (AT) functionals of type 1 and 2, respectively.The corresponding crack profiles are given by dAT- Some authors like [28] prefer w(d) = d 2 because it has a local minimizer at d = 0. Thus, in the absence of mechanical strain, the unfractured solution d ≡ 0 is a minimizer of the total energy.As a result, no additional constraints need to be applied to ensure that d ≥ 0. However, this argument becomes void when solver technology is available that can handle the explicit constraints 0 ≤ d ≤ 1.In contrast, for the AT-1 functional we have w = 0 in the intact state d = 0. Together with the constraint d ∈ [0, 1] this leads to a threshold, i.e., a minimum load required to cause damage [40].A numerical comparison of the AT-1 and AT-2 functionals can be found in [10].
Kuhn et al. [28] proposed to regard the Ambrosio-Tortorelli functionals as special instances of the general family defined by ( 12) . The Ambrosio-Tortorelli functionals are obtained by setting β = 0 for AT-1 and β = −1 for AT-2.Further choices of w are proposed in [40], which also contains a detailed stability analysis for one-dimensional problems.We note the following properties of the functional w in (12): Lemma 2.12.The function w given in (12) has the following properties: (1) It fulfills w(0) = 0 and w(1) = 1.
For the rest of the paper we will assume the w(•) takes the form (12) with β ≤ 0 such that w(•) is guaranteed to be convex and quadratic.

Discretization and the algebraic increment potential
We use a fully implicit discretization in time, and Lagrange finite elements for discretization in space.By using a fully implicit time discretization we retain the variational structure of the problem.Most of this section is spent investigating the properties of the increment functional.
3.1.Time discretization.It is shown in [37] that there is a natural time discretization for (1) that consists of sequences of minimization problems.To simplify the presentation we will not derive the time discretization from an energetic formulation of the time-dependent problem (cf.Remark 2.1) but first discretize the variational formulation (4) and then reformulate the time-discrete problem as a sequence of minimization problems.Let the time interval [0, T ] be subdivided by time points t n , n = 0, 1, 2, . . .and denote by (u n , d n ) ∈ H u0 × H d0 the discrete approximation of (u(t n ), d(t n )).Approximating the time derivative ḋ(t n+1 ) by the backward difference quotient (d n+1 − d n )/τ n for the time step size τ n := t n+1 − t n , and inserting this into the time-continuous variational inequality (4) we obtain the time-discrete variational inequality for (u n+1 , d n+1 ) , using the fact that R is 1-homogeneous, and relabeling (v, ê) to (v, e) yields ( 13) . This is the first-order optimality system for the minimization problem (14) (u n+1 , d n+1 ) := arg min with the increment potential Although the variational inequality (13) is not equivalent to the minimization problem ( 14) because E is not convex, we will use the minimization formulation in the following.
Note that the time step size does not appear in this functional, which means that the model is rate-independent.Note also that the increment potential depends on the previous time step only through the indicator functional.Lemma 3.1.Assume that Γ D,u is non-trivial in the sense that its m−1-dimensional Hausdorff-measure is positive.Then the functional Π τ n+1 is coercive on H 1 u0 × H 1 d0 .Proof.Using the uniform coercivity (P 5 ) of ψ(•, d), w(1) > 0, and w(d for some constant C > 0. Using Korn's inequality for u, the Poincaré inequality for d, and the fact that P ext (t n+1 , u) grows at most linearly we get for another constant C > 0 Proof.Since weak lower semicontinuity of the other terms in Π τ n+1 follows from convexity and lower semicontinuity of the integrands, we only need to consider the non-convex term To this end we note that (15) can be written as J(u, d, ∇u) for Since F is a Carathéodory function, non-negative (and thus uniformly bounded from below), and convex in ξ for all (x, (u, d)) ∈ Ω × (R m × R), it satisfies the assumptions of Theorem 3.4 in [12].Now let (u ν , d ν ) (u, d) be a weakly convergent sequence in H 1 u0 × H 1 d0 .Then, by the compactness of the embedding into As a direct consequence of coercivity and weak lower semicontinuity we get existence of a minimizer of the increment functional: Theorem 3.3.There is a solution to the minimization problem (14), i.e., there exists a global minimizer (14) of the previous section is posed on the pair of spaces H 1 u0 for the displacements and H 1 d0 for the damage variable.Let G be a conforming finite element grid for Ω.We discretize the function spaces by standard first-order Lagrangian finite elements.In order to derive an algebraic form of the discretized increment functional we make use of the standard scalar nodal basis {θ i } M i=1 associated to the grid nodes {p 1 , . . ., p M } =: N ⊂ Ω. Identifying the R m -valued and scalar finite element functions u and d with their coefficient vectors u ∈ R M,m and d ∈ R M , respectively, we write where u i,j = u j (p i ) and For the integration we use three kinds of quadrature rules: Integrals of smooth nonlinear terms over a grid element e are approximated using a higher-order quadrature rule e,h dV , while the integral over the nonsmooth term I [d n ,1] (d) is approximated using the grid nodes p i as quadrature points, which is often referred to as lumping.All polynomial terms are integrated exactly using appropriate quadrature rules.Using these approximations we obtain the algebraic increment functional J := Π τ,G n+1 given by ( 16) .
Here the quadrature rule Ω,h (•) dV is given by with positive weights ω e,α on each element e.Notice that we do not need quadrature weights in the last term of J , because the indicator function only takes values in {0, ∞}.
To elucidate the algebraic structure of J we introduce the linear operator L : Then the first part J 0 of the functional can be written as .
Note that for the price of a more complex index notation, the linear operator L can also be written as a sparse matrix with suitable blocking structure.In this case L(•, •) e,α : (R M,m × R M ) → (S m × R) corresponds to the (e, α)-th sparse row of this matrix.
As an approximation of the boundary conditions from H 1 u0 × H 1 d0 we will consider J on the affine subspace H alg = H alg u0 × H alg d0 where The associated homogeneous subspace of H alg is denoted by H alg 0 .In the following we make the assumption that N ∩Γ D,u contains the vertices of at least one boundary grid face, which ensures a discrete Korn inequality such that ε(u) 0 ≥ C u 1 holds for all u ∈ H alg u0 .Furthermore we introduce the discrete feasible set that additionally incorporates the pointwise irreversibility constraints.

3.3.
Properties of the discrete incremental potential.The convergence property of the TNNMG algorithm heavily rely on the algebraic structure of the problem.Hence we now collect the essential structural properties of the algebraic increment functional J .While stronger properties hold true for some splittings of ψ, we only note the necessary properties shared by all of the proposed splittings.In order to preserve the significant properties in the presence of numerical quadrature, we assume that the quadrature rule e,h f dV can at least integrate the isotropic energy f = |ε(u)| 2 F exactly for any finite element function u.Lemma 3.4.The functional J 0 = A + B has the following properties: To see uniform strong convexity, we note that uniform strong convexity of ψ(•, d) implies that there is some η > 0 such that φ(ε, d) Using the exactness assumption on the quadrature rule we get Since this is a weighted sum of convex functions φ(•, d(q e,α )) with positive weights ω e,α , it is itself convex with respect to u.Thus A(u, d) is the sum of a convex function and the function ε(u) 2 0 .Since the latter is strongly convex on H alg u0 independently of d, the same applies to A(u, d) and (A + B)(•, d).
Finally we note that convexity of g and γ imply convexity of (A + B)(u, •).
The TNNMG algorithm is based on a crucial property called block-separability, which states that the nonsmooth part of the objective functional can be written as a sum, such that the sets of independent variables of the addend functionals are disjoint.We note that J = J 0 + ϕ is of the desired form, with a smooth part J 0 = A + B and a block-separable nonsmooth part which can also be written as the indicator functional ϕ(d) = I K alg d (d) of the feasible set K alg d of the n + 1-th time step.Due to the nonsmoothness of ϕ, the smoothness properties of J 0 do obviously not carry over to the full functional J .Furthermore J is in general not convex as a whole.However we still have the following: Lemma 3.5.The functional J is proper, lower semicontinuous, and coercive on H alg .Furthermore it is convex in u and convex in d.
Proof.Being the indicator function of the closed, nonempty, convex set K alg d it is clear that the separable nonsmooth functional ϕ is convex, proper, and lower semicontinuous.Combining this with smoothness of J 0 we get that J is proper and lower semicontinuous.Similarly, convexity in u and d follows from the corresponding properties of J 0 and ϕ.
Using the uniform coercivity (P 5 ) of ψ(•, d), w(1) > 0, and w(d) ≥ 0 (as in the proof of Lemma 3.1) and the exactness assumption on the quadrature rule (as in the proof of Lemma 3.4) we get for some constant C > 0. Now we can proceed as in the proof of Lemma 3.1 to show coercivity of J .

Truncated Nonsmooth Newton Multigrid for brittle fracture
The Truncated Nonsmooth Newton Multigrid method (TNNMG) is designed to solve nonsmooth block-separable minimization problems on Euclidean spaces.In a nutshell, one step of the TNNMG method consists of a nonlinear Gauß-Seideltype smoother and a subsequent inexact Newton-type correction in a constrained subspace.The nonlinear smoother computes local corrections by subsequent (possibly inexact) solving of reduced minimization problems in small subspaces.As the nonlinear smoother is responsible for ensuring convergence, while the Newton corrections accelerate the convergence, the ingredients of the nonlinear smoother have to be selected carefully.
It is a well known result [17] that nonsmooth Gauß-Seidel-type methods can easily get stuck if the subspace decomposition used to construct localized minimization problems is not aligned with the decomposition induced by the block-separable nonsmooth term.In our case, the nonsmooth term ϕ is separable with respect to the decomposition of unknowns induced by the grid vertices.An additional requirement is that the local minimization problems must be uniquely solvable, which is typically ensured by choosing the decomposition such that the local problems are strictly convex.
In view of these requirements we first decompose the space according to the grid vertices and then with respect to the local uand d-degrees of freedom leading to a decomposition Here the m-dimensional subspace V j,u represents the displacement components at the j-th grid vertex, while the one-dimensional subspace V j,d represents the dcomponent at this vertex.All other components are set to zero in these spaces such that the subspace decomposition can be written as a direct sum.For simplicity we use a plain enumeration of these subspaces in alternating order Notice that with this splitting none of the nonsmooth terms ϕ i in (17) couples across different subspaces.Furthermore, by Lemma 3.4 the restriction of J to any affine subspace (u, d) + V i , i = 1, . . ., 2M is convex.
We will now introduce the TNNMG method.For simplicity we first assume that Γ D,u and Γ D,d are empty and that J 0 is C 2 .Let ν ∈ N 0 denote the iteration number.Given a previous iterate U ν = (u, d) ν ∈ R M,m × R M , one iteration of the TNNMG method consists of the following four steps: The algorithm is easily generalized to non-trivial Dirichlet boundary conditions by leaving out all subspaces associated to Dirichlet vertices during the nonlinear smoothing, and by additionally requiring W ν ⊂ H alg 0 for the linear correction subspace.Then, if the initial iterate satisfies the boundary conditions, i.e., if (u d ) 0 ∈ H alg , the method will only iterate within this affine subspace, which preserves the Dirichlet boundary conditions for all iterates.
The canonical choice for the linear correction step ( 21) is a single linear multigrid step, which explains why the overall method is classified as a multigrid method.If a grid hierarchy is available, then a geometric multigrid method is preferable.Otherwise, a suitable constructed algebraic multigrid step for small-strain elasticity problems will work just as well.In the case of the phase-field brittle-fracture increment functional considered here, the nonlinear presmoothing step takes a large part of the run-time of a single iteration.The convergence speed can therefore be improved considerably by doing a small fixed number (larger than 1) of multigrid steps, without appreciably increasing the time per iteration.Section 4.1 will discuss convergence of the method based on an abstract convergence theory.The abstract theory will be used as a guideline for the discussion of nonlinear smoothers in Sections 4.2.Finally 4.3 will discuss the linear correction in more detail.4.1.Convergence results.The TNNMG method was originally introduced for convex problems where global convergence to global minimizers can be shown [18,19,21].These classical results cannot be applied here, due to the non-convexity of J .As a generalization of previous results, [20] introduced an abstract convergence theory that also covers non-convex problems.In the following we will summarize some results from this work.These will later be used as a guideline for specifying how to solve the local subproblems (20) and the linear correction problem (21).In order to simplify the presentation some of the terminology and notation used in [20] is avoided in favor of a more specific notation adjusted to the algorithm as introduced above.
Theorem 4.1.Let J : R L → R ∪ {∞} be coercive, proper, lower semicontinuous, and continuous on its domain, and assume that J (V + (•)) has a unique global minimizer in V i for all i and each V ∈ dom J .Assume that the inexact local corrections W i are given by W i = M i (W i−1 ) for local correction operators having the properties: (1) Monotonicity: Furthermore assume that the initial iterate is feasible, i.e., U 0 ∈ dom J , and that the linear correction is monotone, i.e., J (U ν+1 ) ≤ J (U ν+ 1  2 ).Then any accumulation point U of (U ν ) is stationary in the sense that Proof.This is Theorem 4.1 in [20].
Remark 4.2.Note that the statement requires both global lower semicontinuity of J and continuity on its domain, because neither property implies the other.In the proof given in [20], lower semicontinuity is needed to show that accumulation points of the iteration have finite energy and are thus feasible, while continuity of J on its domain is needed to show that they are stationary.Now we discuss the application of this theorem to the phase-field brittle-fracture problem.First we interpret the stationarity result.Proposition 4.3.Let J be given by (16) and the subspaces V i by (18) and (19).Then any stationary point U in the sense of ( 22) is first-order optimal in the sense of Proof.The stationarity (22) implies the variational inequalities for each subspace V i .Now let W ∈ dom J .Since the splitting ( 18) is direct one can splitting W uniquely into Using the product structure of dom J we find that Summing up the variational inequalities for those W i we obtain the assertion.
Next we investigate the assumptions of the theorem.First we note that J as given in ( 16) is coercive, proper, and lower semicontinuous by Lemma 3.5.Furthermore, J 0 is continuous and the indicator function ϕ is continuous on its domain.Hence the latter is also true for J = J 0 + ϕ.Subspaces V i with odd index i only vary in u(p i ) such that existence of a unique minimizer of J (W + (•))| Vi follows from the strong convexity of J (•, d) shown in Lemma 3.4.For even i these subspaces are associated to nodal damage degrees of freedom d(p i ).Although J (u, •) is in general only convex, but not strictly convex, the restriction J (W + (•))| Vi to a single node is a strictly convex quadratic functional, which again implies existence of a unique minimizer.Finally the monotonicity J (U ν+1 ) ≤ J (U ν+ 1 2 ) of the linear correction is a direct consequence of the damped update.
It remains to identify proper local correction operators M i satisfying the above assumptions.As a first result we show that solving the local minimization problems (20) exactly leads to a convergent algorithm in the sense given above.Lemma 4.4.Let J be given by (16) and the subspaces V i by (18) and (19).Then the exact local solution operators satisfy the assumptions of Theorem 4.1.
Depending on the damaged energy density ψ, solving the restricted problems exactly may not be practical.As a remedy, it is also shown in [20] that inexact minimization is sufficient as long as it guarantees sufficient decrease of the energy.In fact we do not need W i = M i (W i−1 ) exactly but may relax this to J (W i ) ≤ J (M i (W i−1 )) for a suitable continuous M i .However, sufficient decrease is in general hard to check rigorously.In the following we cite one inexact variant from [20] where sufficient descent is guaranteed a priori.Lemma 4.5.Let J be given by (16) and the subspaces V i by ( 18) and (19).For each subspace V i let C i be a symmetric positive definite matrix that satisfies Then the correction operators Proof.This is Lemma 5.8 in [20].

4.2.
Smoothers for brittle fracture problems.The smoother of the TNNMG method performs a sequence of (inexact) minimization problems in low-dimensional subspaces V i .Different approaches are possible here, implementing different compromises between convergence speed, wall-time per iteration, and ease of programming.As there are two types of degrees of freedom, two types of solvers are needed as well.
4.2.1.Subspaces of displacement degrees of freedom.We first consider the subspaces V i = V (i+1)/2,u for odd i spanned by the m displacement degrees of freedom at the vertex p (i+1)/2 .Noting that elements of this subspace only vary in the displacement component, the minimization problem ( 20) is equivalent to (23) arg min Here, the restricted functional L i (v) := J (W i−1 + (v, 0)) takes the form where we have used (u, d) = W i−1 .The precise nature of L i depends on the type of energy splitting used by the model.If the isotropic splittings ( 6) or ( 7) are used, ( 24) is a strictly convex quadratic functional on a vector space, and can be minimized exactly by solving an m × m system of linear equations.
For the anisotropic splittings ( 8) and ( 9), the functional is still strictly convex and once continuously differentiable.The classical Hesse matrix, however, is not guaranteed to exist.However, by Lemma 3.4, the increment functional is semismooth.This suggests various natural choices for local solvers, such as steepestdescent methods or nonsmooth Newton methods [48].When these are used to solve the local problems (20) exactly, global convergence of the overall TNNMG method follows from Theorem 4.1 and Lemma 4.4 above.
However, as mentioned in the previous section, Theorem 4.1 is more general, and also shows convergence for certain types of inexact local solvers (such as the one in Lemma 4.5).Such a setup can make iterations much faster, while keeping the corresponding deterioration of the convergence rate within acceptable limits.Possible approaches are: • One Newton-type step where the ψ -term in the Hessian J 0 of the differentiable part J 0 of J is replaced by the quadratic upper bound (1+k)ψ 0 , that is, the undamaged St. Venant-Kirchhoff energy density (scaled with (1+k)), which is strictly convex and quadratic.By the construction of the splittings in Section 2.1, this term bounds the degraded energy density (5) for any admissible value of d.We call this approach a preconditioned smoother.• One (or another fixed number of) semismooth Newton steps.
• One gradient step with exact line search.For the first variant, global convergence of the TNNMG solver follows from Lemma 4.5.For the other two, the problem of showing convergence is open.Section 5 will show how the first two choices perform in practice.4.2.2.Subspaces spanned by damage degrees of freedom.For subspaces V i = V i/2,d with even i, i.e., subspaces spanned by the damage degrees of freedom at the vertices p i/2 , the minimization problem ( 20) is equivalent to arg min . For all choices of damage functions g described in Section 2 this is a strictly convex quadratic functional on a closed interval, whose minimizer can be computed directly by computing the unconstrained minimizer and projecting onto the admissible interval.We therefore always assume that these problems are solved exactly.4.3.Linear multigrid corrections.For the linear correction step (21) we need to compute a constrained Newton-type correction at least inexactly.This requires to determine the subspace W ν , to compute the constrained first-and second-order derivatives J (U ν+ 1 2 )| Wν and J (U ν+ 1 2 )| Wν ×Wν on this subspace, and finally to solve the system inexactly.
It is easy to see that the largest subspace W ν where J is differentiable in a neighborhood of U ν+ 1 2 = (u ν+ 1 2 , d ν+ 1 2 ) is given by In this subspace the nonsmooth indicator functional ϕ is identical to zero such that we only need to compute first-and second-order derivatives of the smooth part J 0 , which are then restricted to the degrees of freedom that are allowed to be nonzero in W ν .This can easily be achieved by setting rows and columns of J and J to zero for degrees of freedom not contained in W ν .For all splittings where the degraded density ψ is not C 2 , it is at least locally Lipschitz and semismooth.In this case a generalized second-order derivative J 0 (U ν+ 1 2 ) can be used as a replacement of the classical Hesse matrix making (21) a semismooth Newton step.Such a generalized Hessian can be obtained by the following procedure: The density ψ is piecewise C 2 , and every point (ε, d) where it is not C 2 is at the boundary of several subdomains on which ψ is C 2 .Whenever the second derivative ψ needs to be computed at such a point, one uses instead the second derivative from any of these adjacent subdomains.
To practically compute the first and second derivatives of the degraded elastic energy with the splitting (9) based on eigenvalues, recall that the positive and negative parts ψ ± 0 of that splitting are spectral functions where Eig : S m → R m is the ordered eigenvalue function, and the ψ ± 0 : R m → R are invariant under permutations of their arguments.Such functions are called spectral, and they are once or twice differentiable if and only if the functions ψ ± 0 are once or twice differentiable, respectively.The first and second derivatives of ψ ± 0 can be expressed in closed form in terms of the derivatives of ψ ± 0 .This is shown, for example, in [31,Lemma 3.1] for first-order derivatives, and in [31,Theorem 3.3] for secondorder ones.The expressions involve computing the eigenvector decomposition of the argument of ψ ± 0 .For the inexact solution (25) one step of a classical linear multigrid method can be used.Here we only need to take care that the linear smoother can deal with the non-trivial kernel resulting from constraining the linearization.For a linear Gauß-Seidel smoother this amounts to omitting corrections for rows with zero diagonal entry.When using the TNNMG method for problems with a nonquadratic smooth part such as the phase-field brittle-fracture problem considered here, the smoother is relatively expensive.One can then improve the overall convergence rate by doing more than one multigrid iteration, without increasing the time per iteration.

Numerical examples
In this last section we demonstrate the speed and robustness of the TNNMG solver with two numerical examples.For both of them we study four instances from the family of models discussed in this manuscript: the isotropic and the spectral splittings (( 6) and ( 9), respectively), combined with the AT-1 and AT-2 crack density functionals (w(d) = d and w(d) = d 2 , respectively).
In the experiments we will consider two variants of the TNNMG method with two different nonlinear smoothers, both based on the splitting (19) which alternates displacement and damage degrees of freedom.For the first variant-denoted TNNMG-EX in the following-the smoother will solve the local displacement problems (23) inexactly in the sense that one damped (generalized) Newton step is applied.In the second variant-denoted TNNMG-PRE-the smoother will solve approximate local displacement problems exactly.To this end it employs a single Newton-like step where the ψ -term in the Hessian J 0 of the differentiable part J 0 of J is replaced by the quadratic upper bound (1 + k)ψ 0 .
The constrained quadratic local damage problems are solved exactly by both smoothers.By Lemma 4.5, the TNNMG-PRE smoother satisfies the assumptions of the convergence Theorem 4.1, but the TNNMG-EX smoother does so only for the degraded elastic energy density without a splitting (otherwise its local solution operator is not continuous).For the linear correction step (21) we use three standard V (3, 3) linear geometric multigrid steps with a block Gauß-Seidel smoother operating on the canonical (d + 1) × (d + 1) blocks.The coarse linear problems of the multigrid iteration are solved using the UMFPack sparse direct solver [13].Damage degrees of freedom are truncated when they are less than 10 −10 away from their lower bound.The TNNMG iteration is set to run until the relative degraded energy norm of the correction drops below 10 −7 .The degraded energy norm for the coupled problem combines the energy norm of the degraded linear elasticity problem with the energy norm of the AT-2 crack surface energy density. 4e measure iteration numbers, wall-time, and memory consumption.The TNNMG algorithm and the operator-splitting algorithm used for comparison are implemented in C++ using the Dune libraries5 [5,43].

5.
1.An operator splitting method for phase-field brittle fracture models.We compare the performance of the TNNMG method to an operator-splitting method from the literature.This method is described here in some detail to allow readers to reproduce the results.We choose an operator-splitting method because such methods are widely used, and reported to be robust [2,54].
Starting from an initial displacement u 0 and damage field d 0 the operatorsplitting method alternately repeats the following steps: The iteration is terminated under the same condition as the TNNMG method above.
If the degraded elastic energy ψ(ε, d) is of the type (6), which does not distinguish between compressive and tensile strains, then ( 26) is a linear problem with a symmetric positive definite matrix.Our implementation solves these problems using the CHOLMOD direct sparse solver [11]. 6f the degraded elastic energy incorporates a compressive-tensile split such as ( 9), then J is not a quadratic functional in u.In this case, as suggested by Miehe et al. [35,Equation (66)], we perform a single Newton step at (u ν , d ν ), viz.
where ∇ u J und ∇ 2 uu J are the first and second derivatives of J , respectively, with respect to u only.Because of the splitting, the Hesse matrix ∇ 2 uu J does not depend continuously on u.At points where an entry of ∇ 2 uu J jumps, we simply pick one of the one-sided limits, which are elements of the generalized Jacobian of ∇ u J in the sense of Clarke (cf.Remark 2.10).
For the AT-1 and AT-2 models, the damage subproblem ( 27) is a quadratic minimization problem subject to lower bound constraints The Hesse matrix is For the AT-2 model this Hesse matrix is always positive definite.For the AT-1 model, it is only positive definite if the tensile elastic energy ψ + 0 is positive everywhere, and positive semidefinite otherwise.This did not lead to any problems in the numerical tests done for this article.If necessary, a small amount of L 2 -regularization can be added to increase the robustness.
Following a suggestion by Burke et al. [10], we use a projected Newton method to solve the constrained damage problems.We use the variant proposed by Bertsekas [6], because it is well described and, for the quadratic problems considered here, converges locally quadratically [6, Proposition 4].The method combines projected gradient steps for degrees of freedom in the vicinity of the constraints with Newtontype steps for the other degrees of freedom, and an Armijo-style line search.Consult the original article [6] of Bertsekas for a detailed description.In the notation of that article, we let M be the identity matrix and D k (k being the iteration number) the truncated Hesse matrix.This matrix results from the true Hesse matrix ∇ 2 dd J (u ν+1 , d) by replacing the rows and columns for which the degrees of freedom are close to or at the constraint by the corresponding rows and columns of the identity matrix. 7The maximum truncation distance is ε = 10 −5 .For the line search parameters we set β = 0.5 and σ = 0.49.We solve the linear subproblems with the CHOLMOD solver.
Note that the projected Newton method is closely related to the TNNMG method.Indeed, the TNNMG method could also be applied to solve only the damage problem (27) within an operator splitting loop.The main differences are that TNNMG has a presmoother, and that it does not solve the linear correction problems exactly.choice here.However, direct solvers are widely used in practice, and we therefore consider a comparison with such a solver worthwhile. 7Note that this way of truncating degrees of freedom differs from the one employed by the TNNMG method (the construction of the space Wν ) in Section 4.3.
Exploiting symmetry we only simulate on the upper half of the domain (right).

5.2.
Pure tension test of a notched, symmetric specimen.The first numerical example is a two-dimensional, square-shaped notched specimen of size L × L under tension.Due to symmetry we simulate only its upper half.Geometry and boundary conditions are shown in Figure 2. On the top edge a time-dependent normal displacement ū is prescribed, while the horizontal displacement is left free.The bottom edge of the upper half is clamped vertically for all x > L/2, and fixed vertically and horizontally at the single point (L/2, 0), which is where the initial crack tip is.With increasing normal displacement ū, the preexisting crack opens, and the specimen ruptures suddenly, when a limit load is exceeded.By symmetry, the crack energy is accounted for correctly, even though only one half of the crack profile appears in the simulation result.
The simulations are performed with parameters taken from the corresponding experiment in [36], viz.L = 1 mm, Lamé parameters λ = 121 kN/mm 2 and µ = 80 kN/mm 2 , critical energy release rate g c = 2.7 N/mm, and residual stiffness k = 10 −5 .The phase-field regularization parameter is set to l = 0.03125 mm.We apply the loading in 160 steps, and set the displacement load ū at step i to where e 2 is the canonical basis vector point upwards.We start the evolution with no displacement and no damage anywhere.
For the spatial discretization we use three different uniform grids with 256 × 128 (h 1 ), 512 × 256 (h 2 ), and 1024 × 512 (h 3 ) quadrilateral elements, respectively.These were all constructed by uniform refinement of a grid with 32 × 16 elements, and hence the grid hierarchy for the multigrid solver consists of 4, 5, and 6 levels, respectively.To separate the effect of the discretization parameter h from the modeling parameter l, we explicitly study different element sizes h to illustrate the performance of the algorithm, instead of choosing a single element size proportional to the crack width.Relating the grid resolution to the average fracture with l, the average element edge length corresponds to l/8, l/16, and l/32, respectively.We first consider the model with the isotropic splitting (6) of the elastic energy density, where all elastic strains contribute to the degradation of the material.Figure 3 shows the evolution and displacement-force curves, both for the AT-1 and AT-2 functionals.The force plotted here is the total normal force on the top edge of the specimen.
We first compare iteration numbers.At each loading step, the increment problem is solved starting from the solution of the previous time step until the energy norm of the correction normalized by the energy norm of the previous iterate drops below 10 −7 .The upper row of Figure 4 shows the number of iterations for the TNNMG solver with the exact smoother (TNNMG-EX), with logarithmically scaled vertical axis.As can be seen, the iteration numbers remain essentially bounded independently from the grid resolution.The peak shortly before the 150th load step is where the material ruptures.In this situation, the system becomes highly unstable, which results in a higher number of iterations.
After the rupture, the iteration numbers do depend on the grid resolution.This is atypical for a multigrid solver-presumably it is caused by the fact that, given the particular boundary conditions, the completely ruptured specimen is essentially an ill-posed problem.For comparison, the lower row of Figure 4 shows the iteration numbers of the operator splitting method.One can see that for the first two thirds of the loading history, this method needs less iterations than the TNNMG algorithm, and that iteration numbers are independent from the grid resolution.Recall, however, that TNNMG iterations are much cheaper than operator-splitting iterations, because they are essentially a low fixed number multigrid iterations, whereas each operatorsplitting iteration involves solving two global linear systems.In contrast to the TNNMG method, the number of iterations increases with increasing load, and shortly before the peak iteration numbers are higher.
We do not show the iteration numbers of TNNMG with the inexact smoother (TNNMG-PRE), because they coincide with the results for TNNMG-EX.This is not surprising: As the model uses the isotropic splitting of the elastic energy, the local displacement problems are quadratic, and a single Newton step solves them exactly.TNNMG-PRE uses a preconditioner for those quadratic problems, which in this particular situation is very similar to the actual problem.Therefore, preconditioning here has only a limited impact on the smoothing and thus on the speed AT-1 AT-2 AT-1 AT-2 . 2d example, isotropic energy split: Wall-time per degree of freedom per time step, for grid sizes h 1 , h 2 , h 3 of convergence.However, since the preconditioner is independent of the current iterate, the corresponding local matrices can be precomputed.
Wall-time behavior is discussed next.We plot wall-time per time step for the two multigrid variants TNNMG-EX and TNNMG-PRE, and for the operator-splitting method (Figure 5, again with logarithmic vertical axes).The plots show the time normalized by the number of degrees of freedom.
We see that TNNMG is about 2 to 3 times faster than operator-splitting.The time per degree of freedom stays roughly constant for both methods, independent of the grid resolution.Presumably, the superlinear complexity of the direct solver used in the operator-splitting method only shows for larger grids.Table 1 shows the accumulated normalized run-times for the entire load history.It shows that the speed difference for the h 3 grid accumulates to a factor of about 3.5 for the AT-1 model.For the AT-2 model the difference is only a factor of about 2. This is a bit surprising: The factor is larger for the smaller grids, but TNNMG, for the AT-2 model, slows down considerably when going from the h 2 grid to the h 3 grid.This  is caused by higher iteration numbers throughout the load history, as can be seen in Figure 4.The reason for this behavior is unclear.
Figure 5 and Table 1 also show the wall-times for TNNMG-PRE, the TNNMG variant smoothing with an inexact local solver.It can be observed that using that smoother decreases the computation time by roughly further 5% to 30%.This is the effect of precomputing the local preconditioned Hessians in TNNMG-PRE.
Finally, we point out that the AT-1 model is solved more quickly than the AT-2 one, even though the threshold for damage formation makes it more challenging.6) by the spectral one (9).As the example specimen is loaded in tension we expect few differences to the previous experiments, and indeed, the displacement-load curves (in Figure 6) show only minor differences compared to the ones of the isotropic splitting (Figure 3).The purpose of this test is primarily to assess the cost of the two-dimensional eigenvalue decomposition and its derivatives required for the spectral splitting.Figure 7 shows the iteration numbers per time step for the three methods.As the model is not quadratic in the displacement anymore, we now distinguish between the TNNMG-EX and the TNNMG-PRE smoothers.Not surprisingly though, the iteration numbers are virtually identical.The iterations needed by the operator splitting method, shown in the lowest row of Figure 7, have not changed appreciably either.As an exception to this, the TNNMG-PRE method on the h 3 grid needs a much larger number of iterations at the actual rupture, where the problem is very ill-conditioned.While this is only a single load step, it markedly influences the accumulated wall-times.
Figure 8 shows the time needed to solve the increment problems, again normalized by the number of degrees of freedom.One can see that the overall behavior remains unchanged, but that the time needed by TNNMG goes up a bit, in particular for the AT-2 model.As the number of iterations per load step remains largely unchanged, the observed cost increase is caused by the spectral decompositions performed by the smoother.Interestingly, the operator-splitting solver for the AT-2 problem is a bit faster during the rupturing of the specimen than it was for the isotropically split energy.
When looking at the accumulated run-times shown in Table 2 (still normalized by the number of degrees of freedom), we see that the time needed by the operatorsplitting method remains virtually unchanged.The TNNMG methods have gotten slower, though, and have only a small wall-time lead over operator-splitting.As in the case of the isotropic splitting, moving from the h 2 grid to the h 3 grid increases the time per degree of freedom considerably for the AT-2 model.The reason is unclear.
Using the preconditioned smoother instead of the exact one now leads to a significant decrease in the accumulated wall-time, with time reductions between 10% and 40%.This is caused by the fact that the preconditioned smoother computes much fewer eigenvector decompositions.No such speedup can be seen for the AT-1 model on the h 3 grid, where the preconditioned smoother is actually slower than the exact one.As mentioned, this is because the algorithm needs many more iterations at the load step with the rupture in this case.The reason is unknown.AT-1 AT-2 AT-1 AT-2

5.2.3.
Comparison to an interior-point solver from the literature.To allow further assessment of the solver wall-times, we compare them with results published by Wambacq et al. [49,Chapter 4.1], that use an interior-point solver for a very similar example.There, the authors consider the same problem geometry and loading (quantitatively), and similar material parameters.They model the material with an AT-2 functional and an elastic energy using Amor's volumetric-deviatoric splitting [3], see Section 2.1.2.Note that this splitting is cheaper to evaluate than the eigenvalue-based one whose times are given in Table 2.
As we do, the authors simulate a linearly increasing displacement load until a bit after complete rupture, but they use only 60 load steps compared to our 160 ones.Apparently they simulate the entire square, but with 17 421 vertices, their grid is a bit coarser than our h 1 grid (which has 33 153 vertices).They give cumulative simulation times for four types of solvers, of which three are of operator-splitting type (with different local solvers), and one is a monolithic interior-point solver.For the operator-splitting solvers they report wall-times per degree of freedom between 6.9 ms and 27.6 ms.For a better comparison with our simulation with 160 load steps . Boundary value problem of bending test in three dimensions we multiply these times with 160 60 and obtain values between 18.4 ms and 73.5 ms.Likewise, for the monolithic solver they report a cumulative time per degree of freedom of 82.7 ms (for 60 load steps).Scaled to 160 load steps this gives a value of about 220 ms.These times should be compared to the ones of the AT-2 column of Tables 1 and 2. One can see that our implementations are quite a bit faster.
Note, however, that the two example problems are not exactly the same, that termination criteria differ, and that both hardware and implementation differ as well.Also, the grid used by Wambacq et al. [49] is coarser than ours.Therefore, the comparison should not be over-interpreted.On the other hand, the interiorpoint implementation involves sparse matrix factorizations, and it is therefore to be expected that the run-times deteriorate rapidly for increasing grid resolution, in particular for three-dimensional problems.

5.3.
Bending test of a notched bar in three dimensions.The second example uses a three-dimensional object.In three dimensions, stiffness matrices get denser, and hence direct solvers for global linear systems get more expensive and need considerably more memory.In contrast, the TNNMG memory consumption remains linear in the number of unknowns.Also, eigenvalue decompositions are more expensive for 3 × 3 matrices, and we therefore expect larger run-time differences between the isotropic and the spectrally split model, and between the two TNNMG smoothers.
We consider a bending test for a rectangular bar with a triangular notch.In this setting, the decomposition of the elastic energy density plays a crucial role.Under the given loading, parts of the specimen undergo severe compression, and material models that degrade under such compression will therefore show unphysical results.We test the solvers with the isotropic splitting (6) nevertheless, to obtain an idea of the cost of the spectral splitting.
The example setting is again taken from [36].The geometry and the boundary conditions are visualized in Figure 9.The dimensions of the specimen are L x = 8 mm, L y = 2 mm and L z = 1 mm.Width and height of the triangular notch are l 1 = 0.4 mm and l 2 = 0.2 mm, respectively.We use three unstructured hexahedral grids to discretize the domain,with 1 920, 15 360, and 122 880 elements, respectively.In the figures these are denoted by h 1 , h 2 , and h 3 , respectively.The grids were constructed by 2, 3, and 4 steps of uniform refinement of a coarse grid with 30 elements, respectively, and this refinement history is used by the linear geometric multigrid step of the TNNMG algorithm.The grids are graded a bit towards the expected fracture, and have an average edge length of about l/2.3 (for h 1 ), l/4.6 (for h 2 ) and l/9.2 (for h 3 ) there.
As in [36], the Lamé parameters are set to λ = 121 kN/mm 2 and µ = 80 kN/mm 2 .For the further parameters we set g c = 2.7 N/mm, l = 0.2 mm, and k = 10 −5 .The object is loaded with a time-dependent displacement load in downward direction in a strip of width 1.2 mm in the center of the top surface.In that strip, the surface is fixed in z-direction, but free to move in x-direction.The displacement load is set to ūi = −i • 5 • 10 −3 mm • e 3 for the load steps i = 1, . . ., 13.No adaptive load-stepping as in [36] is necessary, because time discretization and solvers are stable enough for this range of load step sizes.The object is fully clamped at a strip of width 0.4 mm at the left end of the lower surface.At the right end of the same surface, a strip of width 0.4 mm is fixed in the y and z directions.
For each time step, we solve the increment problem with the same solver settings as for the two-dimensional example: Each iteration starts at the solution of the previous load step, and we terminate when the degraded energy norm of the correction, scaled by the degraded energy norm of the previous iterate, drops below 10 −7 .The linear correction step (21) consists of three geometric multigrid V (3, 3)cycle steps with a block-Gauß-Seidel smoother, and damage degrees of freedom are truncated when their current value is less than 10 −10 away from the lower obstacle.5.3.1.Isotropic splitting.We first consider the model with the isotropic splitting (6) of the elastic energy density, i.e., the model where all elastic strains contribute to the degradation of the material.We expect unphysical results: Virtually all damage will happen in the vicinity of the load, whereas the region around the notch will remain intact.Indeed, the simulation results in Figure 10 show that this is exactly what is happening.We are nevertheless interested in the solver behavior for this model, primarily to assess the cost of the spectral splitting in the next section.No costly eigenvalue decompositions are necessary for the isotropic splitting considered here, and the local displacement problems (24) are quadratic.Consequently, the two smoother variants TNNMG-EX and TNNMG-PRE solve almost the same problems, and we expect them to behave more or less the same, too, with a small run-time advantage for TNNMG-PRE.
To assess solver performance, we again first compare iteration numbers.Figure 11 shows the iteration numbers per time step needed by the two TNNMG variants, and by the operator-splitting iteration.Note again that the vertical axis is scaled logarithmically.In contrast to the two-dimensional experiment we did notice different iteration numbers for the exact and preconditioned smoothers, and we therefore show both plots.
We see that TNNMG needs about 10 iterations for each load step for the AT-1 model.For the AT-2 model, for which the damage variable changes throughout the domain immediately, the solver also needs about 10 iterations on the h 2 grid, but almost 10 times as many for the other two grids for the first 5 load steps.Operator-splitting, on the other hand, behaves roughly the same for both models.Unlike TNNMG, it needs less iterations initially, but many more later on.More specifically, the method needs around 10 iterations or even less for the first 5 load steps, but after that it consistently needs about 100 iterations per load step for all grids, with one outlier even going up to 1000 iterations.For an attempt of an explanation, recall that in this example the specimen never breaks into two parts (Figure 10), and the problem therefore remains much better conditioned than the previous ones.In this situation, the monolithic TNNMG method seems to have a clear advantage.
The difference in iteration numbers together with the lower cost of TNNMG iterations compared to operator-splitting iterations leads to a tremendous speed difference.As shown in Figure 12, for all load steps beyond the first few, the TNNMG method needs about two orders of magnitude less wall-time than the operator-splitting method.This difference can also be seen in Table 3, which shows the accumulated wall-times for the entire load history, normalized by the number of degrees of freedom.We see that TNNMG-EX is about 25 to 35 times faster than operator-splitting for the AT-1 model, and 25 to 30 times faster for the AT-2 model.Using the preconditioned smoother makes the wall-time decrease further.In this three-dimensional situation, not having to recompute the matrix block diagonal entries does lead to large savings.The speed advantage of TNNMG-PRE over operator-splitting rises to about 42-58 for the AT-1 model, and to 33-57 for the AT-2 one.
5.3.2.Spectral splitting of the elastic energy density.In the second set of experiments we replace the unsplit degraded energy density ( 6) by the one with the spectral splitting according to (9). Figure 13 shows two snapshots from the problem evolution, and the reaction force as a function of the applied displacement.One can clearly see the differences to the simulation with the unsplit energy.As one would expect, the damage now happens primarily near the notch, and the specimen does now break into two parts.Figure 14 shows the iteration numbers for TNNMG with both smoothers and the operator-splitting method again.Iteration numbers for the multigrid algorithms and both AT models are roughly the same.In particular, TNNMG-PRE, the multigrid method with the inexact smoother avoiding the costly 3 × 3 eigenvalue decomposition, does not need more iterations than TNNMG-EX.The operatorsplitting method needs slightly less iterations than TNNMG for the first few load steps, and a few more after that.
When looking at wall-times, the TNNMG algorithm is consistently faster than the operator-splitting method.Figure 15 shows the normalized times in the same way as in the previous sections.The TNNMG-EX method is roughly five times faster in each load step.This shows the effect of the cheaper iteration times: In three spaces dimensions, tangent matrices are denser than in two dimensions, and direct solving of linear systems with such matrices gets more expensive.Consequently, operator-splitting iterations get relatively more expensive than multigrid ones.This is also reflected in Table 4, which shows the accumulated run-times for the entire load history per degree of freedom.Here, TNNMG-EX is about 4 to 7 times faster than operator-splitting for the AT-1 model, and roughly 5 times faster for the AT-2 model.
The speed difference gets even larger when using the preconditioned smoother.Computing eigenvector decompositions (and the derivative transformation formulas for spectral functions; Section 4.3) is much more expensive for 3 × 3 matrices than for 2 × 2 ones.Therefore, using the preconditioned smoother, which avoids many of these computations, saves a lot of time.As can be seen from Figure 15 and Table 4, TNNMG-PRE is about 7 to 12 times faster for the AT-1 model, and 8 to 10 times faster for the AT-2 model.Finally, we investigate the memory consumption of the two solver implementations.Figure 16 shows the maximum amount of memory used by the two algorithms for the problem of this section, as a function of the grid size.Memory consumption was measured using the valgrind-massif tool 8 with the --pages-as-heap option.One can see that TNNMG clearly requires a linear amount of memory.On the other hand, the proportionality factor is rather large, because TNNMG needs the full Newton tangent matrix, restriction operators, and approximate tangent matrices on coarser grid levels.Also, in order to obtain first-and second-order derivatives efficiently during the local solves, our implementation In contrast, operator splitting does not use coarser grid levels and only assembles tangent matrices for the displacement and damage problems separately.This leads to a smaller memory footprint for small grids.Once the grid gets finer, though, the superlinear memory complexity of the direct solver begins to show.To highlight this effect we measured the memory consumption for one further step of grid refinement.The resulting grid has 983 040 elements, and the simulation using TNNMG needs  .3d example, spectral energy split: Wall-time per degree of freedom over time step, for grid sizes h 1 , h 2 , h 3 about 47 GB of memory.The corresponding operator-splitting implementation, however, would exhaust even a machine with 512 GB of main memory (Figure 16).

( 1 )
D (u,d) E(t, u, d) + ∂ ḋR(d, ḋ) 0, where D (u,d) E(t, u, d) means the Gâteaux derivative with respect to the second and third arguments of E, and ∂ ḋR(d, ḋ) is the convex subdifferential with respect to the second argument of the dissipation potential R.

Figure 1 .
Figure 1.Crack width definition for the AT-2 functional

Figure 3 .
Figure 3. 2d example, isotropic energy split: Evolution at time steps 145 and 160 computed with the TNNMG algorithm, and displacement-force curves for the AT-1 (top) and AT-2 (bottom) models

Figure 6 .
Figure 6.2d example, spectral energy split: Evolution at time steps 145 and 160, and displacement-force curves for the AT-1 (top) and AT-2 (bottom) models

h 2 h 3 (Figure 8 .
Figure 8. 2d example, spectral energy split: Wall-time per degree of freedom per time step, for grid sizes h 1 , h 2 , h 3

Figure 10 .
Figure 10.3d example, isotropic energy split: Configurations at load steps 7 and 13, and displacement-force curves.

Figure 13 .
Figure 13.3d example, spectral energy split: Evolution at time steps 7 and 13, and displacement-force curves precomputes and stores the values and derivatives of the shape functions at all quadrature points (what Miehe et al. call the global interpolation matrix [35, Chap.4.2]).

h 2 h 3 (
Figure 15.3d example, spectral energy split: Wall-time per degree of freedom over time step, for grid sizes h 1 , h 2 , h 3

Figure 16 .
Figure 16.3d example: Memory consumption for solving one increment problem of the AT-2 functional

We remind the readers that the gradient ∇ψ(•, d) is called semismooth if for any point A ∈ S m and any direction V ∈ S m the limit lim n→∞ G n V n exists
and is unique for all sequences V n and G n with V n → V and G

Table 1 .
2d example, isotropic energy split: Total wall time per degree of freedom (in milliseconds) Spectral splitting of the elastic energy density.For the next experiment we exchange the isotropic energy splitting (

Table 3 .
3d example, isotropic energy split: Total wall time per degree of freedom (in milliseconds), for grid sizes h 1 , h 2 , h 3 Figure 11.3d example, isotropic energy split: Iterations per time step, for grid sizes h 1 , h 2 , h 3 Figure 12. 3d example, isotropic energy split: Wall-time per degree of freedom over time step, for grid sizes h 1 , h 2 , h 3

Table 4 .
3d example, spectral energy split: Total wall time per degree of freedom (in milliseconds), for grid sizes h 1 , h 2 , h 3 Figure 14.3d example, spectral energy split: Iterations per time step, for grid sizes h 1 , h 2 , h 3