Computational Mechanics (2023) 72:1059–1089 https://doi.org/10.1007/s00466-023-02330-x ORIG INAL PAPER Truncated nonsmooth Newtonmultigrid for phase-field brittle-fracture problems, with analysis Carsten Gräser1 · Daniel Kienle2 ·Oliver Sander3 Received: 2 December 2022 / Accepted: 24 March 2023 / Published online: 20 May 2023 © The Author(s) 2023 Abstract 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. Keywords Phase-field · Brittle fracture · Spectral strain decomposition · Convex analysis · Nonsmooth multigrid · Global convergence 1 Introduction The equations of phase-fieldmodels of brittle fracture present a number of challenges to the designers of numerical solution algorithms [2]. Even in the small-strain case the equations are nonlinear, due to the multiplicative coupling of the mechani- 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). B Oliver Sander oliver.sander@tu-dresden.de Carsten Gräser graeser@math.fau.de Daniel Kienle kienle@mechbau.uni-stuttgart.de 1 Department Mathematik, Friedrich-Alexander-Universität Erlangen–Nürnberg, Cauerstr. 11, 91058 Erlangen, Germany 2 Institut für Angewandte Mechanik, Universität Stuttgart, Pfaffenwaldring 7, 70569 Stuttgart, Germany 3 Institut für Numerische Mathematik, Technische Universität Dresden, Zellescher Weg 12–14, 01069 Dresden, Germany cal 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 [38] 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 solv- ing these spatial problems is based on operator splitting. Algorithms based on this approach, also known as stag- gered 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 [16], which provides a natural frame- work for convergence proofs. Applications of the operator splitting scheme and its extensions appear, e.g., in [8, 10, 11]. With particular semi-implicit time discretizations it is also possible to solve the spatial problems by solving only one damage problem and one displacement problem [37], without iterating. 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 [16, 17, 52, 54, 55]. For 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s00466-023-02330-x&domain=pdf https://orcid.org/0000-0003-4855-8655 https://orcid.org/0000-0003-1093-6374 1060 Computational Mechanics (2023) 72:1059–1089 the unmodified Newton method only local convergence can be shown, and failure to converge for large load steps is read- ily observed in practice [57]. Therefore, various authors have proposed extensions or modifications of the Newton idea to stabilize the method. In [17] a line search strategy is applied to enlarge the domain of convergence of Newton’s method. [57] and [29] propose to use the BFGS quasi-Newton algo- rithm, claiming that it is more stable than Newton’s method and more efficient than operator splitting. [24] proposed a modified Newton scheme which was later improved by [55] with an adaptive transition from Newton’s method to the modified Newton scheme. Recent results in [26] suggest that nonlinear preconditioning can speed up the solution process tremendously. Finally, the authors of [36, 48] suggest an arc- length method based on the fracture surface, and an adaptive time stepping scheme to enhance the robustness. In sum- mary, while monolithic Newton-type methods are reported to be faster than operator-splitting algorithms, the latter ones are more robust [2, 57]. Various approaches are used in the literature to deal with the damage irreversibility. A natural approach is to regular- ize the constraint, as investigated in [26, 38, 53]. This leads to an additional parameter, and to ill-conditioned tangent matrices [18]. Interior-point solvers implement an auto- matic control of the new parameter; they are investigated in [52]. An alternative approach considers the thermody- namic driving force of the fracture phase-field as a global unknown yielding a three-field formulation which results in a saddle-point principle [38]. A third formulation con- siders 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 [37]. This approach, frequently known as theH-field technique, therefore trades the inequal- ity constraint for the nondifferentiable maximum function. Alternatively, the time discrete problems can be reformulated as semismooth systems by means of so-called comple- mentarity functions. This strategy can be combined with monolithic semismooth Newton techniques [34] or nested Newton and active set methods [24]. Unfortunately, these approaches spoil the variational structure of the spatial prob- lems. Augmented Lagrangian solvers as in [53] introduce extra variables. Closest in spirit to the present manuscript is the use of bound-constrained optimization solvers, used, e.g., in [4, 11, 16, 56]. None of these approaches are fully satisfactory. The effect of the nondifferentiable terms caused by anisotropic splittings of the mechanical energy density as in [4, 38, 49] is rarely discussed in the literature. Hybrid for- mulations 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 New- ton or operator-splitting iteration. This is expensive, even if efficient multigrid methods are used for the linear sub- problems (as, e.g., in [25]). When the methods use direct sparse solvers for the linear tangent problems, memory con- sumption 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 [39]. As a conse- quence, 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 [38] 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 meth- ods are able to solve variational nonsmooth problems from mechanics efficiently without the need for solving global linear systems of equations. In works like [20, 22, 23, 27, 28, 47], such multigrid methods have shown to be vastly more efficient than operator-splitting or Newton-basedmeth- ods. 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 NewtonMultigrid method (TNNMG) can treat the pointwise constraints of the increment problems directly, i.e., without artificial regularization or tricks like the H-field technique [37]. 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 uncon- strained 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 [12]. We cover elastic ener- gieswith various types of anisotropic splittings, including the splittingbasedon strain eigenvalues of [38]. For theproofswe use results from the convex analysis of spectral functions [5, 123 Computational Mechanics (2023) 72:1059–1089 1061 44]. Extension to the slightly more general damage models of [35] is straightforward, provided the stored elastic energy has the required convexity and smoothness properties. In con- trast to the multilevel trust region method proposed in [27], 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: Sect. 2 discusses a framework of small-strain phase-field models for brittle frac- ture, and shows the range of applicability of the TNNMG solver. Section3 introduces the natural fully implicit time discretization, and proves existence of solutions for the spa- tial problems. In both sections we pay particular attention to the mathematical properties of the energy functionals. In Sect. 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 Sect. 5. Our reference for comparison, briefly revisited in Sect. 5.1, is the operator-splitting iteration proposed in [12], which uses a projected Newton method [7] 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. 2 Phase-fieldmodels of brittle fracture This section presents a range of phase-field models for brittle fracture, and discusses its smoothness and convexity proper- ties. 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+∇uT ). Following [38], 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 thiswe select twonot necessarily equal subsets�D,u, �D,d ⊂ ∂� of the domain boundary, and require u = u0 on �D,u, d = d0 on �D,d , where u0 and d0 are two given functions. Displacement and damage field evolve together, gov- erned by a system of coupled nonsmooth partial differential equations. Disregarding inertia effects, we obtain a rate- independent system in the sense of [39]. Such a system can be written using the Biot equation D(u,d)E(t, u, d) + ∂ḋR(d, ḋ) � 0, (1) where D(u,d)E(t, u, d) means the Gâteaux derivative with respect to the second and third arguments ofE , and ∂ḋR(d, ḋ) is the convex subdifferential with respect to the second argu- ment of the dissipation potential R. In this equation, E is a potential energy, which we assume to be of the form E(t, u, d) = ∫ � ψ(ε(u), d) dV + ∫ � gcγ (d,∇d) dV + Pext(t, u) + ∫ � I[0,1](d) dV . (2) The term ψ is a degraded elastic energy density, and will be discussed in detail in Sect. 2.1. The term γ models the local crack surface density, and will be discussed in Sect. 2.2. The number gc is Griffith’s critical energy release rate, a material parameter. Pext represents time-dependent volume and surface forces, which drive the evolution.We assume that Pext is linear and H1(�)-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 IK : R → R ∪ {∞}, IK(x):= { 0 if x ∈ K, ∞ otherwise. For a closed, convex, nonempty set K, the functional IK is convex, lower semicontinuous, and proper. Adding the con- straint 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 practi- cally 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 H1(�,Rm) × H1(�,R). Incorporating the boundary con- ditions leads to the affine subspace H1 u0 × H1 d0 := { v ∈ H1(�,Rm) ∣∣ v|�D,u = u0 } × { v ∈ H1(�) ∣∣ v|�D,d = d0 } . The second term of the Biot equation (1) is ∂ḋR(d, ḋ), where R is the dissipation potential R(d, ḋ) = ∫ � I[0,∞)(ḋ) dV . (3) 123 1062 Computational Mechanics (2023) 72:1059–1089 It implements the pointwise non-healing condition ḋ ≥ 0 as proposed by [38]. Note that R(d, ·) : H1(�) → [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 func- tional (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 equiva- lent to the variational inequality 〈 D(u,d)E(t, u, d), (v, e) − (u, ḋ) 〉 + R(e) ≥ R(ḋ) ∀(v, e) ∈ H1 u0 × H1 0 . (4) Likewise, it is equivalent to the coupled system 〈DuE(t, u, d), v〉 = 0 ∀v ∈ H1 0,〈 DdE(t, u, d), e − ḋ 〉 + R(e) ≥ R(ḋ), ∀e ∈ H1 0 , where we have denoted byH1 0 × H1 0 :=H1 u0 × H1 d0 − (u0, d0) the homogeneous space corresponding to H1 u0 × H1 d0 . Since these two notions of solutions are based on first-order deriva- tives of E andR, they may lead to local minimizers or saddle points of the functional E . To overcome this, there are alter- native 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 [39, 50], where they are discussed for damage-related and more general rate-independent pro- cesses. Remark 2.2 In the engineering literature, the same problem is sometimes formulated as ∂(u̇,ḋ)�(u̇, ḋ; u, d) � 0, with the rate potential �(u̇, ḋ; u, d):= d dt E0(u, d) + Pext(t, u̇) + R(d, ḋ), (see, e.g., [38, Section 4.2]), where E0(u, d) denotes the parts of the energy E(t, u, d) that do not explicitly depend on t . If Pext(t, u) is linear in u, and if the problem is sufficiently smooth, then this formulation is equivalent to the Biot equa- tion (1). Indeed, we then have �(u̇, ḋ; u, d) := D(u,d)E0(u, d)(u̇, ḋ) + Pext(t, u̇) +R(d, ḋ), and ∂(u̇,ḋ)�(u̇, ḋ; u, d) = D(u,d)E0(u, d) + Pext(t, ·) + ∂ḋR(d, ḋ) = D(u,d) [E0(u, d) + Pext(t, u) ] + ∂ḋR(d, ḋ). Requiring this to contain 0 is the Biot equation (1). 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 ψ0(ε) = λ 2 tr[ε]2 + μ tr[ε2], withLamé parametersμ > 0 andλ > − 2 3μ.With this choice of parameters the quadratic functional ψ0 is strongly convex on Sm , the set of real symmetric m × m matrices. The undamaged energy density is split as ψ0(ε) = ψ+ 0 (ε) + ψ− 0 (ε) 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 g : [0, 1] → [0, 1], and the energy density ψ : Sm × [0, 1] → R takes the form ψ(ε, d) = [g(d) + k]ψ+ 0 (ε) + (1 + k)ψ− 0 (ε). (5) The residual stiffness k > 0 guarantees awell-posed problem in case of fracture. Various different degradation functions have appeared in the literature [12, 30, 49]. While the details vary, there appears to be agreement on the following properties: Assumption 2.3 The degradation function g : [0, 1] → [0, 1] is differentiable, monotone decreasing, and fulfills g(0) = 1 and g(1) = 0. 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., [54]). However, in the typ- ical 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. 123 Computational Mechanics (2023) 72:1059–1089 1063 Note that several authors assume g′(1) = 0 in order to ensure that the evolution does not lead to values of d larger than 1 (cf., e.g., [42]). We do not need this assumption on g′ here, because the pointwise constraint d ≤ 1 is enforced explicitly by the energy term (2). The following specific degradation functions all fulfill Assumption 2.3: ga(d) = (1 − d)2 (from [9]) gb(d) = (1 − d)2 · (2d + 1) (from [30]) gc(d) = (1 − d)3 · (3d + 1) (from [30]) gd(d) = exp(bd) − (b(d − 1) + 1) exp(b) (b − 1) exp(b) + 1 , b > 0 (from [49]). Note that the functions ga and gd are strictly convex, but gb and gc 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 litera- ture. 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 fol- lowing essential properties: (P1) ψ(ε, ·) ∈ C2 for all ε ∈ S m andψ(·, d) ∈ C1,1 for all d ∈ [0, 1], i.e., ψ(·, d) is differentiable with locally Lipschitz continuous derivative. (P2) The gradient ∇ψ(·, d) is semismooth for all d ∈ [0, 1]. (P3) The gradient ∇ψ(·, d) is globally Lipschitz con- tinuous uniformly in d, i.e., there exists L ≥ 0 independent of d such that for all matrices A, B ∈ S m we have |∇ψ(A, d) − ∇ψ(B, d)| ≤ L|A − B|F , where | · |F denotes the Frobenius norm.2 (P4) ψ(·, d) : Sm → R is strongly convex uniformly in d, i.e., there exists η > 0 independent of d such that for all matrices A, B ∈ S m we have ψ ( t A + (1 − t)B, d ) ≤ tψ(A, d) + (1 − t)ψ(B, d) −1 2 ηt(1 − t)|A − B|2F . 1 The stress-based splitting of [49] is left for future work. 2 While we only require L ≥ 0 in (P3), the strong convexity (P4) in fact implies L > 0. (P5) ψ(·, d) is coercive uniformly in d in the sense that there exists C > 0 independent of d such that ψ(ε, d) ≥ C |ε|2F .3 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→∞GnVn exists and is unique for all sequences Vn and Gn with Vn → V and Gn ∈ ∂(∇ψ(·, d))(A + tnVn) for a sequence tn ↘ 0. The set ∂(∇ψ(·, d))(A) denotes Clarke’s generalized Jaco- bian of the locally Lipschitz continuous map ∇ψ(·, d) : S m → S m at A ∈ S m (cf. [45]). Notice that the strong convexity (P4) implies strong monotonicity of ∇ψ , i.e., 〈∇ψ(A) − ∇ψ(B), A − B〉 ≥ η|A − B|2F . Here and in the following we denote by 〈·, ·〉 : V∗ ×V → R the duality pairing of a vector space V . For the splittings considered in the following we will only prove (P1) and (P2) directly, and show that the simplified assumptions of the following lemma hold true. This then implies (P3), (P4), and (P5). Lemma 2.5 Let ψ+ 0 and ψ− 0 be convex, non-negative, and differentiable with Lipschitz continuous gradients ∇ψ+ 0 and ∇ψ− 0 . Then ψ satisfies (P3), (P4), and (P5). Proof Let L+ and L− be the Lipschitz constants of∇ψ+ 0 and ∇ψ− 0 , respectively. Then ∇ψ(·, d) is Lipschitz continuous with uniform Lipschitz constant (1+ k)(L+ + L−), because g(d) + k ≤ 1 + k. To show strong convexity, we first note that ψ0 = ψ+ 0 + ψ− 0 is strongly convex on S m with a modulus η > 0 inde- pendent of d. Now consider the function ε �→ψ(ε, d) − kψ0(ε) = g(d)ψ+ 0 (ε) + ψ− 0 (ε). 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 ψ(ε, d) ≥ kψ0(ε) ≥ k η 2 |ε|2F . �� Despite those strong properties of ψ(·, d) we note that ψ(ε, d) is not convex in d and ε together for any of the splittings considered below. 3 Under the additional assumption that 0 = ψ(0, d) ≤ ψ(ε, d) holds for all ε and d one can show that (P4) implies (P5). 123 1064 Computational Mechanics (2023) 72:1059–1089 2.1.1 Isotropic splitting In this model, any strain will lead to damage. The splitting is therefore ψ+ 0 (ε) = ψ0(ε), ψ− 0 (ε) = 0. (6) Without proof, we note the following simple properties of the energy density ψ defined by (5) and this splitting: Lemma 2.6 The energy density ψ defined in (5) with the isotropic splitting (6) has the properties (P1)–(P5). Further- more ψ(·, d) has the stronger property that it is in C∞ and quadratic for all d ∈ [0, 1]. 2.1.2 Volumetric decompositions The isotropic model is of only limited use, because it pro- duces fracturing for all kinds of strain. In [31], Lancioni and Royer-Carfagn obtained better results by letting only the deviatoric strain contribute to the degradation. They intro- duced the split ψ+ 0 (ε) = ψ0(dev ε), ψ− 0 (ε) = ψ0(vol ε), with the deviatoric–volumetric strain splitting vol ε:= tr ε m I , dev ε:=ε − vol ε. With these definitions, the energies are ψ+ 0 (ε) = (μ m + λ 2 ) (tr ε)2, ψ− 0 (ε) = μ ( ε2 − 1 m (tr ε)2 ) = μ dev ε : dev ε. (7) Lemma 2.7 The energy density ψ defined in (5) with the isotropic volumetric splitting (7) has the properties (P1)– (P5). Furthermore ψ(·, d) has the stronger property that it is in C∞ and quadratic for all d ∈ [0, 1]. Proof C∞-smoothness and thus (P1) and (P2) are straight- forward. The fact that ψ+ 0 and ψ− 0 are quadratic, convex, and non-negative allows to derive (P3), (P4), and (P5) from Lemma 2.5 and implies that ψ(·, d) is also quadratic. �� The decomposition of Lancioni andRoyer-Carfagni is still isotropic. [4] 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 x2 = 〈x〉2+ + 〈x〉2−, they proposed the energy split ψ+ 0 (ε) = (μ m + λ 2 ) 〈tr ε〉2+, ψ− 0 (ε) = (μ m + λ 2 ) 〈tr ε〉2− + μ dev ε : dev ε, (8) where only the tensile volumetric strain contributes to dam- age. Lemma 2.8 The energy density ψ defined in (5) with the anisotropic volumetric splitting (8) has the properties (P1)– (P5). Furthermore ψ(·, d) is not C2, unless g(d) = 1. Proof We first note that the squared ramp functions 〈·〉2± are convex, C1,1 with derivatives having a global Lipschitz con- stant 2, and piecewise C2 (in the sense of [51, Definition 2.19]). Hence the functions ψ± 0 are also C1,1 with glob- ally Lipschitz gradients and piecewiseC2, which shows (P1) and (using Lemma 2.5) (P3). Being piecewise C2 implies semismoothness (P2) of ∇ψ(·, d) [51, Proposition 2.26]. Noting that μ/m + λ/2 > 0, convexity of the squared ramp functions furthermore implies that the functions ψ± 0 are also convex and non-negative, which by Lemma 2.5 provides (P4) and (P5). For g(d) + k = 1 the functional ψ(·, d) is quadratic and thus C2. In the case g(d) �= 1, if ψ(·, d) would be C2, then the function t �→ ψ(t I , d) would also be C2. However, this function takes the form ψ(t I , d) = (μ m + λ 2 ) m2t2 { g(d) + k if t ≥ 0, 1 + k if t < 0 and is thus piecewise quadratic but not C2 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 [38]. To define this splitting it is convenient to introduce the ordered eigenvalue function Eig : Sm → R m on the spaceSm of symmetricm×mmatrices,mapping any symmetricmatrix 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 ψ± 0 (ε):=λ 2 〈 m∑ i=1 Eig(ε)i 〉2 ± + μ m∑ i=1 〈Eig(ε)i 〉2±. (9) Note that this indeed defines a splitting ψ0 = ψ+ 0 +ψ− 0 . For this splitting we will make the additional assumption that λ ≥ 0. 123 Computational Mechanics (2023) 72:1059–1089 1065 To quantify the properties of ψ(ε, d) with respect to the strain tensor ε we use the theory of spectral functions [32]. To this end we note that we can write ψ± 0 as ψ± 0 = ψ̂± 0 ◦ Eig : Sm → R with ψ̂± 0 (λ):=λ 2 〈 m∑ i=1 λi 〉2 ± + μ m∑ i=1 〈λi 〉2±. The functions ψ̂± 0 are symmetric in the sense that ψ̂± 0 (λ) does not depend on the order of the entries of λ ∈ R m . Having this formwe 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 (P1)– (P5). Furthermore ψ(·, d) is not C2, unless g(d) = 1. Proof We will first show (P1)–(P5). An essential ingredi- ent is that the squared ramp functions 〈·〉2± are non-negative, piecewise quadratic, and convex. (P1) The squared ramp functions 〈·〉2± and thus ψ̂± 0 are C1,1. Now [44, Proposition 4.3] shows that the spectral func- tions ψ± 0 = ψ̂± 0 ◦ Eig are also C1,1. Hence the same applies to ψ(·, d). (P2) The squared ramp functions 〈·〉2± are piecewise C2 functions. Hence the gradients ∇ψ̂± 0 are piecewise C1 functions (in the sense of [51, Definition 2.19]) and thus semismooth [51, Proposition 2.26]. Now [44, Proposition 4.5] provides semismoothness of∇ψ± 0 and thus of∇ψ(·, d). (P3) Since the functions ψ̂± 0 are piecewise quadratic and C1,1 the gradients ∇ψ̂± 0 are globally Lipschitz continuous. NowCorollary 43 of [5] 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 Sm is equipped with the Frobenius norm. Using Lemma 2.5 this implies uniform Lipschitz continuity of ψ(·, d). (P4),(P5) Since the functions ψ̂± 0 are weighted sums of convex, non-negative squared ramp functions with non- negative weights, they are convex and non-negative them- selves. Convexity of the functions ψ± 0 then follows from [5, Theorem41] while non-negativity of those functions is triv- ial. Now Lemma 2.5 provides (P4) and (P5). To characterize second order differentiability of ψ(·, d) we first consider g(d) = 1. Then ψ(·, d) coincides with the quadratic function (1 + k)ψ0 = (1 + k)ψ+ 0 + (1 + k)ψ− 0 and is thus C2. In the case g(d) �= 1, if ψ(·, d) would be C2, then the function t �→ ψ(t E, d) for the fixed matrix E with Ei j = δ1iδ1 j would also be C2. However, this function takes the form ψ(t E, d) = (λ 2 + 1 ) t2 { g(d) + k if t ≥ 0, 1 + k if t < 0, and is thus piecewise quadratic but not C2 in t = 0. �� Remark 2.10 One can show that Sm decomposes into finitely many disjoint subsets Ai such that ψ(·, d) is twice contin- uously differentiable in the interior of each of these sets. A matrix ε ∈ S m is in the intersection of several Ai 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 Ai . Semismoothness essentially means that such generalized derivatives provide an approximation that can be exploited in a generalized Newton method. Remark 2.11 The additional assumption λ ≥ 0 is essential for convexity of ψ(·, d). To see this we consider for m = 2 the line segment { D(t) = diag(−1, t) ∣∣ t ∈ (0, 1) } ⊂ S 2. Then, along this line segment, ψ(·, 1) is quadratic and takes the form ψ(D(t), 1) = kμt2 + (1 + k) (λ 2 (t − 1)2 + μ ) = ( kμ + (1 + k) λ 2 ) t2 + (1 + k) ( −λt + λ 2 + μ ) 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 [43] γ (d,∇d):= 1 4cw (w(d) l + l|∇d|2 ) , with parameters cw and l, and a parameter function w : [0, 1] → [0, 1]. Motivated by the seminal work of Modica andMortola [40, 41], the associated crack surface functional C(d):= ∫ � γ (d,∇d) dx is called a Modica–Mortola functional. The internal length scale parameter l controls the size of the diffusive zone 123 1066 Computational Mechanics (2023) 72:1059–1089 between a completely intact and a completely damagedmate- rial. For l → 0 the regularized crack surface yields a sharp crack topology in the sense of �-convergence [40, 41]. For a given function w, the normalization constant cw 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 potentialsw, 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 cw:= ∫ 1 0 √ w(s) ds. (10) 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 ∈ � : x1 = 0}, and interpret it as a crack in �. To approximate � by a damage field d : � → [0, 1], assume that d minimizes the crack surface energy C(d):= ∫ � γ (d,∇d) dx, subject to the crack conditions d(0, ξ) = 1 and d(±∞, ξ) = 0 for all ξ ∈ �. We find that d is given by d(x) = d̂(|x1|/l) where d̂ : [0,∞] → [0, 1] solves the normalized scalar Euler–Lagrange equation w′(d̂) − 2d̂ ′′ = 0 in (0,∞), with boundary condition d̂(0) = 1, or, equivalently, the ini- tial value problem, d̂ ′(s) = − √ w(d̂(s)) s > 0, d̂(0) = 1. (11) Using (11) we find that the crack surface energy of the min- imizer d is C(d) = 2|�| 4cw ∫ ∞ 0 w(d̂) + |d̂ ′|2ds = −2|�| 4cw ∫ ∞ 0 2 √ w(d̂)d̂ ′ ds = |�| cw ∫ 1 0 √ w(s) ds. Thus cw 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 pro- file, rescaled by the length parameter l. In order to relate the length parameter to the actual crack width, we use the cone construction from [37] and define the crack width to be the average width of a tangential finite cone fitted to the crack profile and the zero function (Fig. 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̂—theconehaswidth 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 w(d) = d, cw = 2 3 and w(d) = d2, cw = 1 2 . They are referred to in the literature as Ambrosio–Tortorelli (AT) functionals of type 1 and 2, respectively. The corre- sponding crack profiles are given by d̂AT-1(s) = {( 1 − s 2 )2 if s < 2, 0 otherwise and d̂AT-2(s) = exp(−s). Some authors like [30] prefer w(d) = d2 because it has a local minimizer at d = 0. Thus, in the absence of mechani- cal 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 thresh- old, i.e., a minimum load required to cause damage [43]. A numerical comparison of the AT-1 and AT-2 functionals can be found in [12]. Kuhn et al. [30] proposed to regard the Ambrosio– Tortorelli functionals as special instances of the general family defined by w(d) = (1 + β(1 − d))d, (12) with β ∈ [−1, 1]. The Ambrosio–Tortorelli functionals are obtainedby settingβ = 0 forAT-1 andβ = −1 forAT-2. Fur- 123 Computational Mechanics (2023) 72:1059–1089 1067 Fig. 1 Crack width definition for the AT functionals ther choices of w are proposed in [43], 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. (2) It is strictly monotone increasing on [0, 1] for all β ∈ [−1, 1]. (3) It is convex for all β ≤ 0, and strictly convex for all β < 0. 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. 3 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 discretizationwe 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 [39] that there is a natural time discretization for (1) that consists of sequences of minimization prob- lems. 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 tn , n = 0, 1, 2, . . . and denote by (un, dn) ∈ Hu0 × Hd0 the discrete approximation of (u(tn), d(tn)). Approximating the time derivative ḋ(tn+1) by the backward difference quotient (dn+1 − dn)/τn for the time step size τn :=tn+1 − tn , and inserting this into the time-continuous variational inequal- ity (4) we obtain the time-discrete variational inequality for (un+1, dn+1) 〈 D(un+1,dn+1)E(tn+1, un+1, dn+1), (v, e) − ( un+1, (dn+1 − dn)/τn )〉 +R(e) − R( (dn+1 − dn)/τn ) ≥ 0 ∀(v, e) ∈ H1 u0 × H1 0 . Testing with (v, e) = ( un+1 + 1 τn (v̂ − un+1), 1 τn (ê − dn) ) for (v̂, ê) ∈ H1 u0 × H1 d0 , using the fact that R is 1- homogeneous, and relabeling (v̂, ê) to (v, e) yields 〈 D(un+1,dn+1)E(tn+1, un+1, dn+1), (v, e) − (un+1, dn+1) 〉 +R(e − dn) − R(dn+1 − dn) ≥ 0 ∀(v, e) ∈ H1 u0 × H1 d0 . (13) This is the first-order optimality system for the minimization problem (un+1, dn+1):= argmin (ũ,d̃)∈H1 u0 ×H1 d0 �τ n+1(ũ, d̃), (14) with the increment potential �τ n+1(u, d):=E(tn+1, u, d) + R(d − dn) = ∫ � ψ(ε(u), d) dV + ∫ � gcγ (d,∇d) dV + Pext(tn+1, u) + ∫ � I[dn ,1](d) dV . 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. 123 1068 Computational Mechanics (2023) 72:1059–1089 Note that the time step size does not appear in this func- tional, 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. Next we investigate the properties—notably the semicontinuity—of the increment potential �τ n+1. While related results forAmbrosio–Tortorelli-type functionals have already been shown in the seminal work [3], we give proofs here for exactly the range of functionals considered in Sect. 2. 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 H1 u0 × H1 d0 . Proof Using the uniform coercivity (P5) of ψ(·, d), w(1) > 0, and w(d) ≥ 0 for d ∈ [0, 1] we get ∫ � ψ(ε(u), d) dV + ∫ � gcγ (d,∇d) dV ≥ C ∫ � |ε(u)|2F + |∇d|2 dV for some constant C > 0. Using Korn’s inequality for u, the Poincaré inequality for d, and the fact that Pext(tn+1, u) grows at most linearly we get for another constant C > 0 �τ n+1(u, d) ≥ C ( ‖u‖21 + ‖d‖21 − 1 − ( ∫ � d dV )2) + ∫ � I[dn ,1](d) dV ≥ C ( ‖u‖21 + ‖d‖21 − 1 − |�|2 ) , where we have used that the constraint d ∈ [0, 1] implies | ∫ � d dV | ≤ |�| in the second inequality. �� Lemma 3.2 The functional �τ n+1 is weakly lower semicon- tinuous on H1 u0 × H1 d0 . 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 ∫ � ψ(ε(u), d) + I[dn ,1](d) dV . (15) To this end we note that (15) can be written as J (u, d,∇u) for J (u, d, ξ):= ∫ � F ( x, (u(x), d(x)), ξ ) dV and the density F : � × (Rm × R) × R m×m → R ∪ {∞} given by F ( x, (u, d), ξ ) = ψ( 12 (ξ + ξ T ), d) + I[dn ,1](d). Since F is a Carathéodory function, non-negative (and thus uniformly bounded from below), and convex in ξ for all (x, (u, d)) ∈ � × (Rm × R), it satisfies the assumptions of Theorem 3.4 in [14]. Now let (uν, dν)⇀(u, d)be aweakly convergent sequence in H1 u0 × H1 d0 . Then, by the compactness of the embedding into L2(�,Rm × R) we get (uν, dν) → (u, d) in L2(�,Rm × R). Furthermore, the H1(�,Rm × R)-weak convergence of (uν, dν) implies L2(�,Rm×m)-weak convergence of ∇uν ξν :=∇uν⇀∇u=:ξ in L2(�,Rm×m), because (u, d) �→ η(∇u) is in H1(�,Rm × R)′ for each η ∈ L2(�,Rm×m)′. Now Theorem 3.4 of [14] provides lim inf ν→∞ J (uν, dν, ξν) ≥ J (u, d, ξ) = J (u, d,∇u). �� As a direct consequence of coercivity and weak lower semicontinuity we get existence of a minimizer of the incre- ment functional: Theorem 3.3 There is a solution to theminimization problem (14), i.e., there exists a global minimizer (un+1, dn+1) ∈ H1 u0 × H1 d0 of �τ n+1. 3.2 Finite element discretization The increment problem (14) of the previous section is posed on the pair of spaces H1 u0 for the displacements and H1 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 }Mi=1 associated to the grid nodes {p1, . . . , pM }=:N ⊂ �. Identifying the R m-valued and scalar finite element functions u and d with their coefficient vectors u ∈ R M,m andd ∈ R M , respectively, we write u j = M∑ i=1 ui, jθi , d = M∑ i=1 diθi , where ui, j = u j (pi ) and di = d(pi ). 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[dn ,1](d) is approximated using the grid nodes pi as quadrature points, which is often referred to 123 Computational Mechanics (2023) 72:1059–1089 1069 as lumping. All polynomial terms are integrated exactly using appropriate quadrature rules. Using these approximations we obtain the algebraic increment functionalJ :=� τ,G n+1 given by J (u, d) := ∫ �,h (ψ(ε(u), d)) dV + ∫ � gcγ (d, ∇d) dV + Pext(tn+1, u) ︸ ︷︷ ︸ =:J0(u,d) + M∑ i=1 I[dn (pi ),1](di ) ︸ ︷︷ ︸ =:ϕ(d) . (16) Here the quadrature rule ∫ �,h(·) dV is given by ∫ �,h f dV := ∑ e∈G ∫ e,h f dV , ∫ e,h f dV := αmax∑ α=1 f (qe,α)ωe,α with positive weights ωe,α on each element e. Notice that we do not need quadrature weights in the last term ofJ , because the indicator function only takes values in {0,∞}. To elucidate the algebraic structure of J we introduce the linear operator L : (RM,m ×R M ) → ((Sm ×R)αmax)G with L(u, d)e,α:=((ε(u))(qe,α), d(qe,α)) α = 1, . . . , αmax, e ∈ G. Then the first part J0 of the functional can be written as J0(u, d) = ∑ e∈G αmax∑ α=1 ψ(L(u, d)e,α)ωe,α ︸ ︷︷ ︸ =:A(u,d) + ∫ � gcγ (d,∇d) dV + Pext(tn+1, u) ︸ ︷︷ ︸ =:B(u,d) . 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,α : (RM,m × R M ) → (Sm ×R) corresponds to the (e, α)-th sparse row of this matrix. As an approximation of the boundary conditions from H1 u0 × H1 d0 we will consider J on the affine subspace H alg = H alg u0 × H alg d0 where H alg u0 :={ u ∈ R M,m | u(p) = u0(p) ∀p ∈ N ∩ �D,u } , H alg d0 :={ d ∈ R M | d(p) = d0(p) ∀p ∈ N ∩ �D,d } . 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 Kalg:=H alg u0 × ( H alg d0 ∩ Kalg d ) , Kalg d := M∏ i=1 [dn(pi ), 1] that additionally incorporates the pointwise irreversibility constraints. 3.3 Properties of the discrete incremental potential The convergence property of the TNNMG algorithm heav- ily 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 prop- erties shared by all of the proposed splittings. In order to preserve the significant properties in the presence of numeri- cal quadrature, we assume that the quadrature rule ∫ e,h f dV can at least integrate the isotropic energy f = |ε(u)|2F exactly for any finite element function u. Lemma 3.4 The functional J0 = A + B has the following properties: (1) J0(·, d) ∈ C1,1 and J0(u, ·) ∈ C2 for any d ∈ Kalg d and u ∈ R M,m. (2) The gradient ∇J0(·, d) is semismooth for any d ∈ Kalg d . (3) The gradient ∇J0(·, d) is globally Lipschitz continuous uniformly in d. (4) J0(·, d) is strongly convex uniformly in d on Halg u0 . (5) J0(u, ·) is convex on Kalg d for any u ∈ R M,m. Proof The smoothness properties, uniform global Lipschitz continuity, and convexity follow from the corresponding properties of ψ and γ , and from linearity of L. To see uniform strong convexity, we note that uniform strong convexity of ψ(·, d) implies that there is some η > 0 such that φ(ε, d):=ψ(ε, d) − η 2 |ε|2F is convex. Using the exactness assumption on the quadrature rule we get ∑ e∈G αmax∑ α=1 φ(L(u, d)e,α)ωe,α = A(u, d) − η 2 ‖ε(u)‖20. Since this is aweighted sumof convex functionsφ(·, d(qe,α)) 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)‖20. Since the latter is strongly convex on H alg u0 independently of d, the same applies to A(u, d) and (A + B)(·, d). 123 1070 Computational Mechanics (2023) 72:1059–1089 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 function- als are disjoint. We note that J = J0 + ϕ is of the desired form, with a smooth part J0 = A+ B and a block-separable nonsmooth part ϕ(d):= M∑ i=1 ϕi (di ), ϕi (ξ):=I[dn(pi ),1](ξ), (17) which can also be written as the indicator functional ϕ(d) = IKalg d (d) of the feasible set Kalg d of the n + 1-th time step. Due to the nonsmoothness ofϕ, the smoothness properties of J0 do obviously not carry over to the full functional J . FurthermoreJ is in general not convex as a whole. However we still have the following: Lemma 3.5 The functional J is proper, lower semicontinu- ous, and coercive on Halg. Furthermore it is convex in u and convex in d. Proof Being the indicator function of the closed, nonempty, convex set Kalg d it is clear that the separable nonsmooth functional ϕ is convex, proper, and lower semicontinuous. Combining this with smoothness of J0 we get that J is proper and lower semicontinuous. Similarly, convexity in u and d follows from the corresponding properties of J0 and ϕ. Using the uniform coercivity (P5) of ψ(·, d), w(1) > 0, and w(d) ≥ 0 (as in the proof of Lemma 3.1) and the exact- ness assumption on the quadrature rule (as in the proof of Lemma 3.4) we get A(u, d) + B(u, d) − Pext(tn+1, u) ≥ C ∫ � |ε(u)|2F + |∇d|2 dV for some constant C > 0. Now we can proceed as in the proof of Lemma 3.1 to show coercivity of J . �� 4 Truncated nonsmooth Newtonmultigrid 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 nonlin- ear Gauß–Seidel-type smoother and a subsequent inexact Newton-type correction in a constrained subspace. The non- linear smoother computes local corrections by subsequent (possibly inexact) solving of reducedminimization problems in small subspaces. As the nonlinear smoother is responsi- ble 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 [19] 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 nons- mooth 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 u- and d-degrees of freedom leading to a decom- position R M,m × R M = (Rm × R)M = M∑ j=1 ( Vj,u + Vj,d ) . (18) Here the m-dimensional subspace Vj,u represents the dis- placement components at the j-th grid vertex, while the one-dimensional subspace Vj,d represents the d-component 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 V2 j−1 = Vj,u, V2 j = Vj,d , j = 1, . . . , M . (19) 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) + Vi , i = 1, . . . , 2M is convex. We will now introduce the TNNMGmethod. For simplic- ity we first assume that �D,u and �D,d are empty and that J0 is C2. Let ν ∈ N0 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: (1) Nonlinear presmoothing (a) Set W0 = Uν (b) For i = 1, . . . , 2M compute W i ∈ W i−1 + Vi as W i ≈ argmin W∈W i−1+Vi J (W) (20) 123 Computational Mechanics (2023) 72:1059–1089 1071 (c) Set Uν+ 1 2 = W2M (2) Inexact linear correction (a) Determine the maximal subspaceWν ⊂ R M,m ×R M such that the restriction J |Wν is C 2 at Uν+ 1 2 (b) Compute cν ∈ Wν as an inexact Newton step on Wν cν ≈ −(J ′′(Uν+ 1 2 )|Wν×Wν )−1(J ′(Uν+ 1 2 )|Wν ) (21) (3) Projection Compute the Euclidean projection cν pr = PdomJ−Uν+1/2 × (cν), i.e., choose cν pr such that U ν+ 1 2 + cν pr is closest to Uν+ 1 2 + cν in domJ (4) Damped update (a) Compute a ρν ∈ [0,∞) such that J (Uν+ 1 2 + ρνcν pr) ≤ J (Uν+ 1 2 ) (b) Set (u, d)ν+1 = Uν+1 = Uν+ 1 2 + ρνcν pr 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 (ud)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 amultigridmethod. 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 func- tional considered here, the nonlinear presmoothing step takes a large part of the run-time of a single iteration. The con- vergence 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 theorywill be used as a guideline for the discussion of nonlinear smoothers in Sect. 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 [20, 21, 23]. These classical results cannot be applied here, due to the non-convexity of J . As a gen- eralization of previous results, [22] 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 cor- rection problem (21). In order to simplify the presentation some of the terminology and notation used in [22] is avoided in favor of a more specific notation adjusted to the algorithm as introduced above. Theorem 4.1 Let J : RL → R ∪ {∞} be coercive, proper, lower semicontinuous, and continuous on its domain, and assume that J (V + (·)) has a unique global minimizer in Vi for all i and each V ∈ domJ . Assume that the inexact local corrections W i are given by W i = Mi (W i−1) for local correction operators Mi : domJ → domJ , Mi − Id : domJ → Vi having the properties: (1) Monotonicity: J (Mi (V )) ≤ J (V ) for all V ∈ domJ . (2) Continuity: J ◦ Mi is continuous. (3) Stability: J (Mi (V )) < J (V ) if J (V ) is not minimal in V + Vi . Furthermore assume that the initial iterate is feasible, i.e., U0 ∈ domJ , and that the linear correction is monotone, i.e., J (Uν+1) ≤ J (Uν+ 1 2 ). Then any accumulation pointU of (Uν) is stationary in the sense that J (U) ≤ J (U + V ) ∀V ∈ Vi , ∀i . (22) Proof This is Theorem 4.1 in [22]. �� 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 [22], lower semicontinuity is needed to show that accumu- lation 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 LetJ be given by (16) and the subspaces Vi by (18) and (19). Then any stationary point U in the sense of (22) is first-order optimal in the sense of 〈∇J0(U),W − U〉 ≤ 0 ∀W ∈ domJ . 123 1072 Computational Mechanics (2023) 72:1059–1089 Proof The stationarity (22) implies the variational inequali- ties 〈∇J0(U),W i − U〉 ≤ 0 ∀W i ∈ domJ ∩ (U + Vi ) for each subspace Vi . Now let W ∈ domJ . Since the split- ting (18) is direct one can splitting W uniquely into W = U + 2M∑ i=1 V i , V i ∈ Vi . Using the product structure of domJ we find that W i := U + V i ∈ domJ ∩ (U + Vi ). Summing up the variational inequalities for those W i we obtain the assertion. �� Next we investigate the assumptions of the theorem. First we note thatJ as given in (16) is coercive, proper, and lower semicontinuous by Lemma 3.5. Furthermore, J0 is continu- ous and the indicator function ϕ is continuous on its domain. Hence the latter is also true for J = J0 + ϕ. Subspaces Vi with odd index i only vary in u(pi ) 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 free- dom d(pi ). 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 mono- tonicity 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 operatorsMi 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 Vi by (18) and (19). Then the exact local solution operators Mi (W):= argmin V∈W+Vi J (V ) satisfy the assumptions of Theorem 4.1. Proof This is Lemma 5.1 in [22]. �� Depending on the damaged energy density ψ , solving the restricted problems exactly may not be practical. As a rem- edy, it is also shown in [22] that inexact minimization is sufficient as long as it guarantees sufficient decrease of the energy. In fact we do not needW i = Mi (W i−1) exactly but may relax this to J (W i ) ≤ J (Mi (W i−1)) for a suitable continuous Mi . However, sufficient decrease is in general hard to check rigorously. In the following we cite one inex- act variant from [22] where sufficient descent is guaranteed a priori. Lemma 4.5 Let J be given by (16) and the subspaces Vi by (18) and (19). For each subspace Vi let Ci be a symmetric positive definite matrix that satisfies 〈∇J0(W + V ) − ∇J0(W), V 〉 ≤ 〈CiV , V 〉 ∀W ∈ domJ , V ∈ Vi . Then the correction operators Mi (W):= argmin V∈(W+Vi )∩domJ J (W) + 〈∇J0(W), V − W〉 +1 2 〈Ci (V − W), V − W〉 satisfy the assumptions of Theorem 4.1. Proof This is Lemma 5.8 in [22]. �� 4.2 Smoothers for brittle fracture problems The smoother of the TNNMG method performs a sequence of (inexact) minimization problems in low-dimensional sub- spaces Vi . Different approaches are possible here, imple- menting 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 Vi = 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 prob- lem (20) is equivalent to argmin (v,0)∈Vi Li (v). (23) Here, the restricted functional Li (v):=J (W i−1 + (v, 0)) takes the form Li (v) = ∫ �,h (g(d) + k)ψ+ 0 (ε(u + v)) + (1 + k)ψ− 0 (ε(u + v)) dV + Pext(tn+1, v) + const, (24) where we have used (u, d) = W i−1. The precise nature of Li 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. 123 Computational Mechanics (2023) 72:1059–1089 1073 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 steepest-descentmethods or nonsmooth Newton methods [51]. When these are used to solve the local problems (20) exactly, global convergence of the overall TNNMGmethod follows from Theorem 4.1 and Lemma 4.4 above. However, as mentioned in the previous section, Theo- rem 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 conver- gence rate within acceptable limits. Possible approaches are: • One Newton-type step where the ψ ′′-term in the Hessian J ′′ 0 of the differentiable part J0 of J is replaced by the quadratic upper bound (1 + k)ψ ′′ 0 , that is, the undam- aged St.Venant–Kirchhoff energy density (scaled with (1 + k)), which is strictly convex and quadratic. By the construction of the splittings in Sect. 2.1, this termbounds 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 prob- lem of showing convergence is open. Section5 will show how the first two choices perform in practice. 4.2.2 Subspaces spanned by damage degrees of freedom For subspaces Vi = Vi/2,d with even i , i.e., subspaces spanned by the damage degrees of freedom at the vertices pi/2, the minimization problem (20) is equivalent to argmin (0,v)∈Vi Li (v), with a restricted functional Li (v):=J (W i−1 + (0, v)). For all choices of damage functions g described in Sect. 2 this is a strictly convex quadratic functional on a closed interval, whoseminimizer 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 cν ≈ −(J ′′(Uν+ 1 2 )|Wν×Wν )−1(J ′(Uν+ 1 2 )|Wν ) (25) at least inexactly. This requires to determine the subspaceWν , to compute the constrainedfirst- 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 Wν = { U = (u, d) ∣∣ di = 0 if (dν+ 1 2 )i /∈ (dn(pi ), 1) } . In this subspace the nonsmooth indicator functionalϕ is iden- tical to zero such that we only need to compute first- and second-order derivatives of the smooth part J0, which are then restricted to the degrees of freedom that are allowed to be nonzero inWν . This can easily be achieved by setting rows and columns ofJ ′ andJ ′′ to zero for degrees of freedom not contained inWν . For all splittingswhere the degraded density ψ is not C2, it is at least locally Lipschitz and semismooth. In this case a generalized second-order derivativeJ ′′ 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 C2, and every point (ε, d) where it is not C2 is at the boundary of several subdomains on which ψ is C2. Whenever the second derivative ψ ′′ needs to be com- puted 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 ψ± 0 = ψ̂± 0 ◦ Eig : Sm → R, where Eig : Sm → R m is the ordered eigenvalue function, and the ψ̂± 0 : Rm → 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 [33, Lemma 3.1] for first-order derivatives, and in [33, Theorem3.3] for second-order ones. The expressions involve computing the eigenvector decomposition of the argument of ψ± 0 . 123 1074 Computational Mechanics (2023) 72:1059–1089 For the inexact solution (25) one step of a classical lin- ear multigrid method can be used. Here we only need to take care that the linear smoother can deal with the non-trivial ker- nel 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 problemswith 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. 5 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) = d2, respectively). In the experiments we will consider two variants of the TNNMG method with two different nonlinear smoothers, both basedon 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 prob- lems exactly. To this end it employs a single Newton-like step where the ψ ′′-term in the Hessian J ′′ 0 of the differen- tiable part J0 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 con- vergence 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 continu- ous). 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 itera- tion are solved using the UMFPack sparse direct solver [15]. Damage degrees of freedom are truncated when they are less than 10−10 away from their lower bound. The TNNMG iter- ation 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.4 We 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 libraries 5 [6, 46]. 5.1 An operator splittingmethod 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, 57]. Starting from an initial displacement u0 and damage field d0 the operator-splitting method alternately repeats the fol- lowing steps: uν+1 = argmin u∈H alg u0 J (u, dν) (26) dν+1 = argmin d∈H alg d0 J (uν+1, d). (27) The iteration is terminated under the same condition as the TNNMG method above. For simplicity we omit any modifi- cations such as proposed, e.g., in [10]. If the degraded elastic energy ψ(ε, d) is of the type (6), which does not distinguish between compressive and ten- sile strains, then (26) is a linear problem with a symmetric positive definite matrix. Our implementation solves these problems using the CHOLMOD direct sparse solver [13].6 If 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 [37, Equation (66)], we perform a single Newton step at (uν, dν), viz. uν+1 = uν − ( ∇2 uuJ (uν, dν) )−1∇uJ (uν, dν), where ∇uJ und ∇2 uuJ are the first and second derivatives of J , respectively, with respect to u only. Because of the splitting, the Hesse matrix ∇2 uuJ does not depend continu- ously on u. At points where an entry of ∇2 uuJ jumps, we 4 Because the AT-1 model does not induce a norm. 5 www.dune-project.org 6 It is well-known that sparse direct solvers do not scale well to larger problems both run-time and memory-wise. Therefore, an alternative method such as linear multigrid may be a better choice here. However, direct solvers are widely used in practice, and we therefore consider a comparison with such a solver worthwhile. 123 www.dune-project.org Computational Mechanics (2023) 72:1059–1089 1075 simply pick one of the one-sided limits, which are elements of the generalized Jacobian of ∇uJ in the sense of Clarke (cf. Remark 2.10). For the AT-1 and AT-2 models, the damage subprob- lem (27) is a quadraticminimization problem subject to lower bound constraints dν i ≤ di ≤ 1 i = 1, . . . , M . The Hesse matrix is ( ∇2 ddJ (uν+1, d) ) i j = ∫ � [ g′′(d)θiθ jψ + 0 (ε(uν+1)) + gc 2cwl θiθ j ︸ ︷︷ ︸ AT-2 only + gcl 2cw ∇θi∇θ j ] dx, with g′′(d) = 2 for our choice g(d) = (1−d)2. 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 L2- regularization can be added to increase the robustness. Following a suggestion by [12], we use a projected New- ton method to solve the constrained damage problems. We use the variant proposed by [7], because it is well described and, for the quadratic problems considered here, converges locally quadratically [7, Proposition 4]. The method com- bines projected gradient steps for degrees of freedom in the vicinity of the constraints with Newton-type steps for the other degrees of freedom, and an Armijo-style line search. Consult the original article [7] of Bertsekas for a detailed description. In the notation of that article, we let M be the identity matrix and Dk (k being the iteration number) the truncated Hesse matrix. This matrix results from the true Hesse matrix ∇2 ddJ (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.7 The 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. 7 Note that this way of truncating degrees of freedom differs from the one employed by the TNNMG method (the construction of the space Wν ) in Sect. 4.3. 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 Fig. 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 dis- placement ū, the preexisting crack opens, and the specimen ruptures suddenly, when a limit load is exceeded. By symme- try, 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 [38], viz. L = 1mm, Lamé parameters λ = 121kN/mm2 and μ = 80kN/mm2, critical energy release rate gc = 2.7N/mm, and residual stiff- ness k = 10−5. The phase-field regularization parameter is set to l = 0.03125mm. We apply the loading in 160 steps, and set the displacement load ū at step i to ūi = i · 2 · 10−5 mm · e2, i = 1, . . . , 160, where e2 is the canonical basis vector pointing upwards. We start the evolution with no displacement and no damage any- where. For the spatial discretization we use three different uni- form grids with 256 × 128 (h1), 512 × 256 (h2), and 1024 × 512 (h3) quadrilateral elements, respectively. These were all constructed by uniform refinement of a grid with 32×16 elements, and hence the grid hierarchy for the multi- grid 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 width l, the average element edge length corresponds to l/8, l/16, and l/32, respectively. 5.2.1 Isotropic splitting We first consider the model with the isotropic splitting (6) of the elastic energy density, where all elastic strains con- tribute to the degradation of the material. Figure3 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. 123 1076 Computational Mechanics (2023) 72:1059–1089 Fig. 2 Notched square with a vertical displacement load (left). Exploiting symmetry we only simulate on the upper half of the domain (right) Fig. 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 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 Fig. 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 rup- tures. 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— 123 Computational Mechanics (2023) 72:1059–1089 1077 Fig. 4 2d example, isotropic energy split: iterations per time step, for grid sizes h1, h2, h3 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 Fig. 4 shows the iter- ation 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 res- olution. Recall, however, that TNNMG iterations are much cheaper than operator-splitting iterations, because they are essentially a low fixed number multigrid iterations, whereas each operator-splitting iteration involves solving two global linear systems. In contrast to the TNNMGmethod, the num- ber 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 sin- gle 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 of convergence. How- ever, since the preconditioner is independent of the current iterate, the corresponding localmatrices can be precomputed. Wall-time behavior is discussed next. We plot wall-time per time step for the twomultigrid variants TNNMG-EX and TNNMG-PRE, and for the operator-splitting method (Fig. 5, againwith 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-splittingmethodonly 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 h3 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-2model, slows down considerably when going from the h2 grid to the h3 grid. This is caused by higher iteration numbers through- out the load history, as can be seen in Fig. 4. The reason for this behavior is unclear. 123 1078 Computational Mechanics (2023) 72:1059–1089 Fig. 5 2d example, isotropic energy split: wall-time per degree of freedom per time step, for grid sizes h1, h2, h3 Table 1 2d example, isotropic energy split: total wall time per degree of freedom (in milliseconds) AT-1 AT-2 TNNMG-EX TNNMG-PRE OS TNNMG-EX TNNMG-PRE OS h1 5.89 5.67 20.47 6.74 5.39 29.19 h2 9.56 6.58 20.51 9.94 7.18 29.43 h3 8.13 7.81 28.84 20.78 18.29 40.92 Figure 5 and Table 1 also show the wall-times for TNNMG-PRE, the TNNMGvariant smoothingwith an inex- act 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 precondi- tioned 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. 5.2.2 Spectral splitting of the elastic energy density For the next experiment we exchange the isotropic energy splitting (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 Fig. 6) show only minor differences compared to the ones of the isotropic splitting (Fig. 3). The purpose of this test is primarily to assess the cost of the two-dimensional eigenvalue decomposition and its deriva- tives required for the spectral splitting. Figure7 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 opera- tor splitting method, shown in the lowest row of Fig. 7, have not changed appreciably either. As an exception to this, the TNNMG-PRE method on the h3 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. 123 Computational Mechanics (2023) 72:1059–1089 1079 Fig. 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 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 iter- ations per load step remains largely unchanged, the observed cost increase is caused by the spectral decompositions per- formed by the smoother. Interestingly, the operator-splitting solver for theAT-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 free- dom), we see that the time needed by the operator-splitting method remains virtually unchanged. The TNNMGmethods 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 h2 grid to the h3 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-1model on the h3 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. 5.2.3 Comparison to an interior-point solver from the literature Toallow further assessment of the solverwall-times,we com- pare them with results published by [52, Chapter 4.1], that use an interior-point solver for a very similar example. There, the authors consider the same problem geometry and load- ing (quantitatively), and similar material parameters. They model the material with an AT-2 functional and an elastic energy using Amor’s volumetric–deviatoric splitting [4], see Sect. 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 dis- placement 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 17421 vertices, their grid is a bit coarser than our h1 grid (which has 33153 ver- tices). They give cumulative simulation times for four types 123 1080 Computational Mechanics (2023) 72:1059–1089 Fig. 7 2d example, spectral energy split: iterations per time step, for grid sizes h1, h2, h3 Table 2 2d example, spectral energy split: total wall time per degree of freedom (in milliseconds) AT-1 AT-2 TNNMG-EX TNNMG-PRE OS TNNMG-EX TNNMG-PRE OS h1 10.20 9.09 20.99 14.78 11.17 29.59 h2 19.58 11.70 21.76 23.55 16.78 29.68 h3 16.26 21.79 32.84 39.64 28.65 40.96 123 Computational Mechanics (2023) 72:1059–1089 1081 Fig. 8 2d example, spectral energy split: wall-time per degree of freedom per time step, for grid sizes h1, h2, h3 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.9ms and 27.6ms. For a better comparison with our simulation with 160 load steps we multiply these times with 160 60 and obtain values between 18.4ms and 73.5ms. Likewise, for the monolithic solver they report a cumulative time per degree of freedom of 82.7ms (for 60 load steps). Scaled to 160 load steps this gives a value of about 220ms. These times should be com- pared 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,8 and that both hardware and implementation differ as well. Also, the grid used by [52] is coarser than ours. Therefore, the compar- ison should not be over-interpreted. On the other hand, the interior-point implementation involves sparse matrix factor- izations, and it is therefore to be expected that the run-times deteriorate rapidly for increasing grid resolution, in particular for three-dimensional problems. 8 The termination criterion in [52] involves Lagrange multipliers for a KKT system that are not present in the problem formulation considered here. 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 getmore expensive and need considerably more memory. In contrast, the TNNMG mem- ory consumption remains linear in the number of unknowns. Also, eigenvalue decompositions aremore 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 [38]. The geome- try and the boundary conditions are visualized in Fig. 9. The dimensions of the specimen are Lx = 8mm, Ly = 2mm and Lz = 1mm.Width and height of the triangular notch are l1 = 0.4mm and l2 = 0.2mm, respectively. We use three 123 1082 Computational Mechanics (2023) 72:1059–1089 Fig. 9 Boundary value problem of bending test in three dimensions unstructured hexahedral grids to discretize the domain, with 1920, 15360, and 122880 elements, respectively. In the fig- ures these are denoted by h1, h2, and h3, respectively. The grids were constructed by 2, 3, and 4 steps of uniform refine- ment of a coarse grid with 30 elements, respectively, and this refinement history is used by the linear geometric multi- grid 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 h1), l/4.6 (for h2) and l/9.2 (for h3) there. As in [38], the Lamé parameters are set to λ = 121kN/mm2 and μ = 80kN/mm2. For the further param- eters we set gc = 2.7N/mm, l = 0.2mm, and k = 10−5. The object is loaded with a time-dependent displacement load in downward direction in a strip of width 1.2mm 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 displace- ment load is set to ūi = −i · 5 · 10−3 mm · e3 for the load steps i = 1, . . . , 13. No adaptive load-stepping as in [38] 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.4mm at the left end of the lower surface. At the right end of the same surface, a strip of width 0.4mm 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 cor- rection, scaled by the degraded energy norm of the previous iterate, drops below10−7. The linear correction step (21) con- sists of three geometric multigrid V (3, 3)-cycle steps with a block-Gauß–Seidel smoother, and damage degrees of free- Fig. 10 3d example, isotropic energy split: configurations at load steps 7 and 13, and displacement–force curves 123 Computational Mechanics (2023) 72:1059–1089 1083 Fig. 11 3d example, isotropic energy split: iterations per time step, for grid sizes h1, h2, h3 Table 3 3d example, isotropic energy split: total wall time per degree of freedom (in milli- seconds), for grid sizes h1, h2, h3 AT-1 AT-2 TNNMG-EX TNNMG-PRE OS TNNMG-EX TNNMG-PRE OS h1 4.15 2.63 109.77 13.31 11.94 393.92 h2 6.11 3.74 216.42 7.12 4.07 223.96 h3 10.34 6.28 357.35 14.36 9.40 349.97 123 1084 Computational Mechanics (2023) 72:1059–1089 Fig. 12 3d example, isotropic energy split: wall-time per degree of freedom over time step, for grid sizes h1, h2, h3 dom 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 elas- tic strains contribute to the degradation of the material. For this particular benchmark we expect unphysical results: Vir- tually all damage will happen in the vicinity of the load, whereas the region around the notch will remain intact. Indeed, the simulation results in Fig. 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 prob- lems (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 iter- ation numbers. Figure11 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 imme- diately, the solver also needs about 10 iterations on the h2 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 manymore later on.More specifically, the method needs around 10 iterations or even less for the first 5 load steps, but after that it consistently 123 Computational Mechanics (2023) 72:1059–1089 1085 Fig. 13 3d example, spectral energy split: evolution at time steps 7 and 13, and displacement–force curves 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 (Fig. 10), and the problem there- fore 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 Fig. 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 num- ber 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. In situations such as hydraulic fracturing where an isotropic model is justified, the advan- tage of TNNMG is therefore considerable. 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 onewith the spectral split- ting according to (9). Figure13 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. Iter- ation 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 operator-splitting 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. Fig- 123 1086 Computational Mechanics (2023) 72:1059–1089 Fig. 14 3d example, spectral energy split: iterations per time step, for grid sizes h1, h2, h3 ure15 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, tan- gent matrices are denser than in two dimensions, and direct solving of linear systemswith suchmatrices getsmore expen- sive. 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 decom- positions (and the derivative transformation formulas for spectral functions; Sect. 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 com- putations, saves a lot of time. As can be seen from Fig. 15 123 Computational Mechanics (2023) 72:1059–1089 1087 Fig. 15 3d example, spectral energy split: wall-time per degree of freedom over time step, for grid sizes h1, h2, h3 Table 4 3d example, spectral energy split: total wall time per degree of freedom (in milliseconds), for grid sizes h1, h2, h3 AT-1 AT-2 TNNMG-EX TNNMG-PRE OS TNNMG-EX TNNMG-PRE OS h1 23.52 14.69 124.63 23.42 12.87 128.93 h2 38.26 21.36 275.54 46.02 25.07 217.47 h3 64.76 36.55 266.59 78.29 38.31 393.45 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. Figure16 shows the maximum amount of memory used by the two algorithms for the prob- lem of this section, as a function of the grid size. Memory consumption was measured using the valgrind-massif tool 9 with the--pages-as-heapoption.One can see that TNNMG clearly requires a linear amount of memory. On the other hand, the proportionality factor is rather large, because 9 https://valgrind.org/ 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 deriva- tives efficiently during the local solves, our implementation precomputes and stores the values and derivatives of the shape functions at all quadrature points (what Miehe et al. call the global interpolation matrix [37, Chap.4.2]). In contrast, operator splitting does not use coarser grid levels and only assembles tangent matrices for the displace- ment and damage problems separately. This leads to a smaller memory footprint for small grids. Once the grid gets finer, 123 https://valgrind.org/ 1088 Computational Mechanics (2023) 72:1059–1089 Fig. 16 3d example: memory consumption for solving one increment problem of the AT-2 functional though, the superlinear memory complexity of the direct solver begins to show. To highlight this effect we mea- sured the memory consumption for one further step of grid refinement. The resulting grid has 983040 elements, and the simulation using TNNMG needs about 47GB of memory. The corresponding operator-splitting implementation, how- ever, would exhaust even a machine with 512GB of main memory (Fig. 16). Funding Open Access funding enabled and organized by Projekt DEAL. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adap- tation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indi- cate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, youwill need to obtain permission directly from the copy- right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. References 1. Alberti G (2000) Variational models for phase transitions, an approach via �-convergence. In: Buttazzo G, Marino A, Murthy MKV (eds) Calculus of variations and partial differential equa- tions: topics on geometrical evolution problems and degree theory, pp 95–114. Springer, Berlin. https://doi.org/10.1007/978-3-642- 57186-2_3 2. Ambati M, Gerasimov T, Lorenzis LD (2015) A review on phase- field models of brittle fracture and a new fast hybrid formulation. Comput Mech 55:383–405. https://doi.org/10.1007/s00466-014- 1109-y 3. Ambrosio L, Tortorelli VM (1990) Approximation of functional depending on jumps by elliptic functional via �-convergence. Commun Pure Appl Math 43(8):999–1036. https://doi.org/10. 1002/cpa.3160430805 4. Amor H, Marigo J-J, Maurini C (2009) Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. JMechPhysSolids 57(8):1209–1229. https://doi.org/ 10.1016/j.jmps.2009.04.011 5. Baes M (2007) Convexity and differentiability properties of spec- tral functions and spectral mappings on Euclidean Jordan algebras. Linear Algebra Appl 422(2):664–700. https://doi.org/10.1016/j. laa.2006.11.025 6. Bastian P, Blatt M, Dedner A, Dreier N-A, Engwer C, Fritze R, Gräser C, Grüninger C, Kempf D, Klöfkorn R, Ohlberger M, Sander O (2021) The DUNE framework: basic concepts and recent developments. Comput Math Appl 81:75–112. https://doi.org/10. 1016/j.camwa.2020.06.007 7. Bertsekas DP (1982) Projected Newton methods for optimiza- tion problems with simple constraints. SIAM J Control Optim 20(2):221–246. https://doi.org/10.1137/0320018 8. Bourdin B (2007) Numerical implementation of the variational formulation for quasi-static brittle fracture. Interfaces Free Bound 9(3):411–430 9. Bourdin B, Francfort G, Marigo JJ (2000) Numerical experiments in revisited brittle fracture. J Mech Phys Solids 48:797–826 10. Brun MK, Wick T, Berre I, Nordbotten JM, Radu FA (2020) An iterative staggered scheme for phase field brittle fracture propa- gation with stabilizing parameters. Comput Methods Appl Mech Eng. 361:11275