M A Badri

G Rastiello
email: giuseppe.rastiello@cea.fr

E Foerster

Preconditioning strategies for vectorial finite element linear systems arising from phase-field models for fracture mechanics

Keywords: Preconditioning, Algebraic multigrid, Phase-field, Brittle fracture, Vectorial finite element method

Phase-field models are frequently adopted to simulate fracture mechanics problems in the context of the finite element method. To depict fracture, this method involves solving a coupled set of Helmholtz-like damage-field equation and augmented linear momentum balance equation. Solutions to these coupled equations are then used as descriptions of crack propagation phenomena within solids. However, this method imposes a constrain of using extremely fine meshing for properly predicting cracks. For practical problems of interest, this very often leads to linear systems with large sizes that have to be repetitively assembled and solved. As such, iterative solution procedures such as the Krylov subspace based methods for solving these large linear systems within the framework of serial/parallel computing environments become mandatory to obtain results in a feasible time. In this work, the vectorial finite discretization for a hybrid phase-field formulation -a monolithic solving scheme -is presented. The underlying nonlinearity present in the coupled set of equations of the the hybrid phase-field model is dealt through Picard iteration that helps to preserve the symmetry of the linearized system to solve. Due to the symmetric positive definite nature of the finite element linear systems obtained for this problem, the conjugate gradient method makes a standard choice of iterative solution algorithm. Within this article, to improve convergence rates, consequently time to solution, of the conjugate gradient method applied to crack propagation problems, different preconditioning strategies are analyzed, tuned, and discussed. Brittle fracture benchmarks are used to measure the performance of preconditioners which are then applied to massively parallel simulations with millions of unknowns. A series of numerical experiments show that the algebraic multigrid preconditioner is well suited for solving the phase-field model for fracture, being superior to the Jacobi and the block Jacobi preconditioning in all regards: ease of solving the problem, iterations to converge, time to solution, and parallel scaling on more than a thousand processes.

Introduction 1

Fracture mechanics problems are often modeled using phase-field models (see e.g., Francfort and Marigo 2 (1998), [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF], [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF], Miehe et al. (2010b)). Mainly due to their strong 3 theoretical backgrounds and robust numerical implementation in the context of the finite element method 4 (FEM), phase-field formulations are, today, regarded among the robust techniques to model fracture me-5 chanics. Particularly when crack initiation, propagation, and branching need to be modeled accurately. Over 6 the past decades, the phase-field formulations have been developed to model diverse problems of fracture mechanics, e.g., to model brittle fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF], Miehe et al., 2010a,b, Ambati et al., 2014), ductile fracture [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity: Variational formulation and main properties[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF], large-deformation fracture (Miehe andSchänzel, 2014, Borden et al., 2018), dynamic fracture (Borden et al., 2012, Hofacker and[START_REF] Hofacker | Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation[END_REF], cohesive fracture (Verhoosel andde Borst, 2013, May et al., 2015), thermomechanical fracture [START_REF] Kuhn | Phase field simulation of thermomechanical fracture[END_REF], to cite but a few. It is backed by experimental evidence [START_REF] Ambati | A phase-field model for ductile fracture at finite strains and its experimental verification[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Pham | Experimental validation of a phase-field model for fracture[END_REF]) that these models can reliably predict the phenomena of fracture. A recent review article (Wu et al., 2018) provides a general overview of different phase-field models available for solving fracture mechanics problems.

Within the framework of phase-field method, the linear momentum balance equation is supplemented by a Helmholtz-like diffusion equation (damage-field equation) describing damage evolution. The resulting system of partial differential equations are discretized and solved simultaneously. In other words, to predict fracture, for a problem undergoing load variations, the scalar damage-field unknown d = d ∈ [0, 1] is solved simultaneously with the vector displacement field u until the internal body forces are in equilibrium with the external applied forces. The resulting field d then serves as an indicator to the state of the solid, i.e., d = 1 indicates cracked state and d = 0 indicates otherwise. In a finite element context, phase-field formulations lead to regularize damage evolution through introducing an internal length scale such that damage is diffused over a band of elements, which is finite in width. In that sense, phase-field models share some features with other non-local formulations (Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Bažant | Nonlocal Integral Formulations of Plasticity and Damage: Survey of Progress[END_REF][START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF][START_REF] Giry | Stress-based nonlocal damage model[END_REF][START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF], and in particular with gradient-damage formulations [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF]Nedjar, 1996, Peerlings et al., 2001). Comparative studies between phase-field and gradient-damage formulations were presented in (de Borst andVerhoosel, 2016, Mandal et al., 2019).

Although the solution procedure seems fairly straightforward, prerequisite of extremely fine meshing for capturing cracks, and large-scale simulations lead to computationally challenging scenarios. For standard numerical methods, like the finite element method (FEM), one thus ends up solving linear systems with millions of unknowns. Further, nonlinear coupling and (pseudo) time-dependent simulations, requires repeated solutions of these linear systems. Such reasons highly motivate the application of preconditioned iterative solution algorithms to solve the arising linear systems in a serial or a parallel computing framework.

The primary aim of the current article is to identify and highlight different preconditioning options to improve convergence, consequently time to solution, for the finite element linear systems that arise in the phase-field modeling strategy for fracture. Besides this, we also propose a monolithic way of approximating the hybrid phase-field model via the vectorial finite element method which provides straightforward means of parallelization for coupled set of equations [START_REF] Badri | Vectorial finite elements for solving the radiative transfer equation[END_REF].

Literature pertaining use of preconditioners in finite element phase-field modeling of fracture mechanics is scarce. Particularly, till date no study on preconditioning strategies for the hybrid phase-field models for fracture has been presented. Nevertheless, preconditioned staggered FEM solution approach for phase-field model was used by [START_REF] Bilgen | A phase-field approach to conchoidal fracture[END_REF] to model conchoidal brittle fracture. The authors investigated the Jacobi and the geometric multigrid preconditioner for the Krylov subspace conjugate gradient (CG) method [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF], in conclusion the geometric multigrid preconditioner outran the Jacobi preconditioner. Much recently, the Krylov subspace GMRES [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]) (for non-symmetric linear systems) preconditioned with geometric multigrid was used for solving brittle fracture problems in the framework of monolithic FEM solution approach [START_REF] Jodlbauer | Matrix-free multigrid solvers for phase-field fracture problems[END_REF]. Since geometric multigrid preconditioning was used in [START_REF] Bilgen | A phase-field approach to conchoidal fracture[END_REF], [START_REF] Jodlbauer | Matrix-free multigrid solvers for phase-field fracture problems[END_REF], the authors ended up using structured meshes for the numerical experiments that appear in these studies. A modified staggered FEM solution approach was proposed by [START_REF] Farrell | Linear and nonlinear solvers for variational phase-field models of brittle fracture[END_REF] to solve brittle fracture problems. The authors ended up using the complete LU solver and the Krylov subspace MINRES method [START_REF] Paige | Solution of sparse indefinite systems of linear equations[END_REF] (for symmetric positive indefinite linear systems) preconditioned with algebraic multigrid to speed up the solver convergence. [START_REF] Heister | Parallel solution, adaptivity, computational convergence, and open-source code of 2D and 3D pressurized phase-field fracture problems[END_REF] solved pressurized phase-field brittle fracture problems using the Krylov subspace GMRES (for non-symmetric linear systems) preconditioned with a block diagonal preconditioner. For solving these block diagonal matrices authors used algebraic multigrid as a standalone solver.

Within this study, the underlying coupled partial differential equations of the hybrid phase-field model by [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF] are solved by utilizing the vectorial FEM. This formulation was chosen for the sake of studying the solving techniques for phase-field models, but different formulations could have been used.

Within the framework of vectorial FEM, the unknown fields u and d, are coupled together in a single vectorial field w = [u, d] and a single vectorial equation then needs to be solved. Two main approaches -staggered approaches [START_REF] Bourdin | Numerical implementation of the variational formulation for quasi-static brittle fracture[END_REF], Miehe et al., 2010a[START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Bilgen | A phase-field approach to conchoidal fracture[END_REF] and monolithic approaches [START_REF] Gerasimov | A line search assisted monolithic approach for phase-field computing of brittle fracture[END_REF][START_REF] Liu | Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model[END_REF][START_REF] Wick | Modified Newton methods for solving fully monolithic phase-field quasi-static brittle fracture propagation[END_REF][START_REF] Jodlbauer | Matrix-free multigrid solvers for phase-field fracture problems[END_REF] -have been used in the literature to handle nonlinearity in the governing equations of the phase-field model. In this study, the Picard iteration method is used to handle the nonlinearity between the coupled partial differential equations of the hybrid phase-field model. Alongside the vectorial FEM this leads to a monolithic approach.

Consequently, the linear systems obtained are symmetric, sparse, and positive definite in nature. Naturally, for memory considerations, the Krylov subspace-based CG method makes a standard choice of iterative solution algorithm for such linear systems. For improving the performance of the CG method, the most effective preconditioner suggested so far is a geometric multigrid, when solving brittle fracture problems on regular structured grids [START_REF] Bilgen | A phase-field approach to conchoidal fracture[END_REF].

Often for many practical problems of interest, unstructured meshes are preferred over structured ones. The reason is that unstructured meshes allow for better geometric representation of complex or simple geometries. Despite unstructured meshing is simpler to perform compared to a structured one, the recommended geometric multigrid preconditioners are not straightforward to set up in that case. Generating a hierarchy of unstructured meshes which are mandatory for setting up the geometric multigrid preconditioners is not a trivial task. It can be argued that the range of applicability of the geometric multigrid is, therefore, limited [START_REF] Xu | Algebraic multigrid methods[END_REF]). An appealing alternative for problems solved with unstructured meshes is another class of multilevel preconditioners, dubbed as algebraic multigrid (AMG). A review by [START_REF] Stüben | A review of algebraic multigrid[END_REF] provides an overview of this method. With the AMG methods, the coarser grids and associated transfer operators are generated algebraically by using the assembled coefficients of linear system itself. One of the motives for this study, is to tune, apply, and report the performance of AMG preconditioning with application to phase-field models and compare it to other common one-level preconditioners such as the Jacobi and the Block Jacobi.

The rest of the paper is organized as follows. In section 2, the boundary value problem of phase-field fracture mechanics model is introduced briefly. Section 3 then presents the vectorial finite element discretization procedure for the phase-field model. The iterative solution technique of CG along with different preconditioning options is presented in section 4. Results and discussion of the time performance, memory usage, convergence behavior, and parameter tuning of different preconditioners is presented next in section 5. The article then ends with some conclusions in section 6, which highlights the main results of this article and provides recommendations on the use of preconditioners for phase-field fracture mechanics modeling.

Mathematical formulations

In this section, we present the mathematical formulation of the hybrid phase-field model [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] for fracture that is used in this work. Such a model is "hybrid" [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF], Wu and Nguyen, 2018) in the sense that it combines features of the so-called "isotropic" [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] and "anisotropic" [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF], Miehe et al., 2010b, Wu et al., 2018) phase-field models. Here, terms isotropic and anisotropic should not be intended according to their meanings in the Continuum Mechanics theory, but following the nomenclature that is generally used in the phase-field community (Miehe et al., 2010a,b). In particular, the term anisotropic designates models such that special splits of the energy/strain/stress are introduced to represent different dissipative processes occurring in the material under prevailing tensile/compressive states.

Problem setting

We consider a n-dimensional solid body Ω ⊂ R n (n = 1, 2, or 3) submitted to quasi-static external loading such that damage occurs. According to the phase-field model, on each point x(x, y, z) ∈ Ω and at any pseudo-time t ∈ [0, T ], the state of the system is defined by two fields: the vector-valued displacement field

u = u(x, t) : Ω × [0, T ] → R n and the scalar-valued damage field d = d(x, t) : Ω × [0, T ] → [0, 1] ⊂ R. From
now on, the dependency on the space and time variables will be omitted for the sake of compactness. Also,

we denote u := [u 1 , u 2 , u 3 ] T
, where u 1 , u 2 , and u 3 correspond to displacements in x, y, and z directions, respectively.

Boundary value problem of hybrid phase-field model

Consider an arbitrary open bounded domain Ω with imposed tractions on Neumann boundary ∂Ω N ⊂ R n , and provided with essential Dirichlet conditions on boundary ∂Ω D ⊂ R n . The split of the domain boundary ∂Ω ⊂ R n is such that ∂Ω = ∂Ω D ∪ ∂Ω N and ∅ = ∂Ω D ∩ ∂Ω N , with over-line denoting a closure. Under quasi-static conditions and in the absence of body/volume forces within Ω, the hybrid phase-field model poses the following boundary value problem (BVP):

find u : Ω × [0, T ] → R n and d : Ω × [0, T ] → [0, 1] such that ∀ t ∈ [0, T ] div (1 -d) 2 σ(u) = 0 ∀ x ∈ Ω, (1)
G c l 0 + 2H + (u) d -G c l 0 ∆d -2H + (u) = 0 ∀ x ∈ Ω, (2)
given the boundary conditions

σ(u, d) • n = t ∀ x ∈ ∂Ω N , (3) u = u ∀ x ∈ ∂Ω D , ( 4
) ∇d • n = 0 ∀ x ∈ ∂Ω N , ( 5
)
where div(•) is the divergence operator applied to (•), ∆(•) = div(∇(•)) denotes the Laplace operator, G c ∈ R + quantifies the material fracture toughness, and l 0 ∈ R + is the length-scale parameter to control the influence of gradient term. n = n(x) : Ω → R n is the unit normal vector for surface ∂Ω, t = t(x, t) :

Ω × [0, T ] → R n denotes imposed tractions, and u = u(x, t) : Ω × [0, T ] → R n imposed displacements. σ(u) = ∂ψ(u)
∂ε denotes the elastic/effective Chauchy's stress tensor with ψ being the elastic energy defined by

ψ(u) = 1 2 λ tr(ε(u))I + µε(u) : ε(u), (6)
here λ ∈ R + and µ ∈ R + are the Lamé parameters, I is the second order identity tensor, tr(•) denotes the trace operator applied to (•), and ε(u) : Ω × [0, T ] → R n×n denotes the second order small strain tensor, i.e., the symmetric part of the gradient of the displacement field ∇u, such that ε(u) = ∇u + ∇ T u /2.

To account for the experimental evidence that cracking is mainly controlled by tensile stress/strain states, [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] proposed to consider that damage diffusion (in eq. ( 2)) be driven by the maximum tensile elastic energy history function

H + (u) : Ω × [0, T ] → R + : H + (u) = max τ ∈[0,t] ψ + (u), (7)
with tensile elastic energy ψ

+ (u) : Ω × [0, T ] → R being defined as ψ + (u) = (1/2)λ tr ε(u) + + µ tr ε 2 + (u)
, where the bracket operator is expanded as

• + = (• + | • |)/2 and ε + (u) =: Ω × [0, T ] → R n×n
denotes the positive part of stain tensor ε. Definition (7) ensures avoiding cracking in compressed regions.

Finally, to prevent crack interpenetration the following constrain is applied

∀x : ψ + (u) < ψ -(u) ⇒ d := 0, (8)
with

ψ -(u) : Ω × [0, T ] → R denoting the compressive elastic energy, ψ -(u) = (1/2)λ tr ε(u) -+ µ tr ε 2 -(u) , where • -= (• -| • |)/2 and ε -(u) : Ω × [0, T ] → R n×n is now the negative part of strain tensor ε 1 .
While the chosen hybrid phase-field model suggests performing elastic energy split (to avoid cracking in compression regions) and makes use of the tensile elastic energy ψ + (in eq. ( 2) via the history function H + ) to control the damage evolution seldom in tensile regions of the domain, there exists other approaches in the literature. Among others, the commonly used ones involve decomposition of elastic energy based on volumetric and deviatoric contributions [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] and spectral decomposition (Miehe et al., 2010a). However, following these approaches one needs to deal with (strong) non-linear system of equations for elasticity (eq. ( 1) becomes nonlinear due to strain tensor decomposition). One of the key advantages of using hybrid phase-field model (and perhaps due to which such model was proposed by [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF])

is that the momentum balance eq. ( 1) is retained in its linear form, thereby reducing its computational complexity. The only nonlinearity to be dealt with in hybrid phase-field formulation is the one due to variational inequality.

Additionally, as the hybrid phase-field model avoids energy/strain/stress splitting for eq. ( 1), the operator

(1 -d) 2 ∂ψ + (u)/∂ε does not exist in the momentum balance eq. ( 1), which is not the case for other phase-field models [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF], Miehe et al., 2010a). It is well known that this operator is non-convex in both variables u and d and leads to difficulties in designing an efficient solution algorithm (Wick, 2017, Kopaničáková and[START_REF] Kopaničáková | A recursive multilevel trust region method with application to fully monolithic phase-field models of brittle fracture[END_REF]. As such, the hybrid phase-field model becomes an interesting option for solving fracture mechanics with ease. The hybrid phase-field model has been proven to work for a wide variety of fracture mechanics problems [START_REF] Doan | Hybrid phase field simulation of dynamic crack propagation in functionally graded glass-filled epoxy[END_REF][START_REF] Jeong | Phase field modeling of crack propagation under combined shear and tensile loading with hybrid formulation[END_REF][START_REF] Hirshikesh | A FEniCS implementation of the phase field method for quasi-static brittle fracture[END_REF]. While this model ensures accurate crack predictions, it is known to reduce the computational cost (due to avoiding strong nonlinearity) by about one order of magnitude compared with other anisotropic phase-field formulations [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF].

Finite element formulations

This section proceeds with the finite element spatial discretization procedure for the BVP discussed in the previous section. A three-dimensional (3D) formulation (n = 3) is considered for the sake of illustration, but identical equations hold for other dimensions.

Discretized variational formulations

On a meshed domain Ω h ∈ Ω ⊂ R n , for eq. ( 1), the mixed finite element variational formulation in the Lagrangian framework for searching the unknown nodal displacements vector

u h = [u 1 , u 2 , u 3 ] T reads: search u h ∈ V h that satisfies ∀ t ∈ [0, T ] : Ω h (1 -d h ) 2 + κ σ(u h ) : ε(v h ) dv = ∂Ω h N t • v h ds ∀ v h ∈ V h , ( 9
)
where κ 1 is a model parameter to prevent numerical singularity when d → 1. In this formulation, the notation ":" is used for the double contraction between tensors (i.e., component-wise tensor product) and V h is a mixed third order vector-valued finite element functional space to approximate vector test function v h and vector trial function u h :

V h = u h ∈ [H 1 (Ω h )] 3 ∀t ∈ [0, T ] | ∀x ∈ ∂Ω h D u h = u , (10)
with H 1 (Ω h ) denoting a square integrable Sobolev functional space. Similarly, for eq. ( 2), the standard finite element variational formulation for the unknown damage scalar d h reads:

search d h ∈ V h that satisfies ∀ t ∈ [0, T ] : Ω h G c l 0 + 2H + (u h ) d h θ h dv + Ω h G c l 0 ∇d h • ∇θ h dv = Ω h 2H + (u h )θ h dv ∀ θ h ∈ V h , ( 11
)
where, V h denotes the scalar finite element functional space to approximate scalar test function θ h and scalar trial function d h :

V h = d h ∈ H 1 (Ω h ) ∀t ∈ [0, T ] d h ∈ [0, 1] . ( 12
)

Vectorial FEM for the hybrid phase-field model

Using appropriate basis functions to approximate u h and d h , eqs. ( 9) and (11) then lead to two separate linear systems A u x u = b u and A d x d = b d . One can consider solving these linear systems repeatedly in a staggered fashion (one after the another). This way of solving the phase-field model has been used several times in the past, e.g., in [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], [START_REF] Bilgen | A phase-field approach to conchoidal fracture[END_REF], [START_REF] Molnár | 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF]Gravouil (2017), Hirshikesh et al. (2018). With the hybrid phase-field model, the staggered solution approach leads to linear equilibrium (with d h frozen) and damage diffusion (with u h frozen) equations to solve. As a consequence, the only iterations needed are those ensuring the coupling between the two equations.

Alternate to this staggered solving method, another elegant is a the monolithic approach. Here, it involves a reformulation of eqs. ( 9) and ( 11) in such a way that a single linear system Ax = b arises and is solved repeatedly to handle nonlinearity in an iterative way. Studies such as Vignollet et al. (2014), Gerasimov andDe Lorenzis (2016), [START_REF] Liu | Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model[END_REF], [START_REF] Kopaničáková | A recursive multilevel trust region method with application to fully monolithic phase-field models of brittle fracture[END_REF] convey that monolithic solution strategy (due to faster convergence) in comparison with the staggered strategy lead to a faster solution time and higher accuracy for a problem. However, it should be pointed out that staggered strategies are more memory efficient than the monolithic ones. As monolithic solvers need to handles a larger matrix system (Ax = b) as opposed to two relatively smaller matrix systems (A u x u = b u and A d x d = b d ) handled by a staggered solver.

Discretized variational formulation

Considerable gains in assembly and solving times can be achieved if a monolithic approach for the hybrid phase-field model is constructed by using the vectorial FEM discretization procedure. For that purpose, we fully couple the finite element fields u h and d h , and solve these by using vectorial finite elements or in other words fully coupled mixed finite element discretization procedure.

To write the vectorial variational formulation for the hybrid phase-field model, we first exploit the arbitrariness and independency of trial fields v h and θ h , and sum up the two variational eqs. ( 9) and ( 11) to obtain:

search (u h , d h ) ∈ V h × V h that satisfies ∀ t ∈ [0, T ] : Ω h 1 + d h 2 + κ σ(u h ) : ε(v h ) dv + Ω h G c l 0 + 2H + (u h ) d h θ h dv + Ω h G c l 0 ∇d h • ∇θ h dv = ∂Ω h N t • v h ds + Ω h 2H + (u h )θ h dv ∀ (v h , θ h ) ∈ V h × V h . ( 13
)
Then, we introduce a fourth order vectorial trial function

w h = [w h 1 , w h 2 , w h 3 , w h 4 ]
T and a fourth order vectorial test function q h = [q h 1 , q h 2 , q h 3 , q h 4 ] T , with w h ∈ W h and q h ∈ W h . Here, W h denotes a fourth order vectorial finite element space:

W h = w h ∈ [H 1 (Ω h )] 4 ∀t ∈ [0, T ] ∀x ∈ ∂Ω h D {w h i } 3 i=1 = ū ∀x ∈ Ω h w h 4 ∈ [0, 1] . ( 14
)
In fact, within the vectorial trial function w h the first three components {w h i } 3 i=1 represent the displacement fields {u h i } 3 i=1 and the fourth component w h 4 represents d h , i.e., the scalar damage field 2 . The vectorial test function q h follows the same logic.

Nonlinear solving based on Picard iterations

The system of equations resulting from (13) is nonlinear. Different approaches are possible to handle this nonlinearity. A usual choice consists in using a Newton-Raphson iterative algorithm. This implies linearizing the variational formulation with respect to the unknown displacement and damage fields, and then solve for their variations between successive iterations. Linearization process leads, however, to a non-symmetric stiffness-like matrix A, which makes CG method non applicable. Here, to obtain a symmetric formulation and preserve a certain similarity with the staggered solution of the chosen hybrid phase-field formulation, choice was made to use Picard iterations (also commonly known as the method of successive substitution).

Using w h and q h , at a given t, the vectorial variational formulations for the hybrid phase-field model can be written down in the following canonical form:

given w h,j ∈ W h search w h,j+1 ∈ W h that satisfies :

Ω h 1 -w h,j 4 2 + κ λ∇({w h,j+1 i } 3 i=1 ) • ∇({q h i } 3 i=1 ) + 2µε({w h,j+1 i } 3 i=1 ) : ε({q h i } 3 i=1 ) dv + Ω h G c l 0 + 2H + ({w h,j i } 3 i=1 ) w h,j+1 4 q h 4 dv + Ω h G c l 0 ∇w h,j+1 4 • ∇q h 4 dv = ∂Ω h N t • [{q i } 3 i=1 ] T ds + Ω h 2H + ({w h,j i } 3 i=1 )q h 4 dv ∀ q h ∈ W h . ( 15
)
In this equation to deal with the nonlinearity, the superscript indices 'j' and 'j +1' are introduced (discretize pseudo-time). Associated to a finite element variable, these indices correspond to the solutions at previous and current iterations of the solving process.

2 For a two-dimensional problem (n = 2) the vectorial trial function w h is of third order,

w h = [w h 1 , w h 2 , w h 3 ]
T where the first two components {w h i } 2 i=1 represent the displacement fields {u h i } 2 i=1 and the third component w h 3 represents the scalar damage field d h .

Finally, to perform a phase-field simulation for all t ∈ [0, T ], at a given t linear systems Ax = b are assembled and solved consecutively. In this study, piecewise linear polynomial basis [P 1 , P 1 , P 1 , P 1 ] T are used to approximate the Sobolev spaces H 1 (Ω h ) 4 in order to assemble the matrix A and the vector b from the left-hand side and right-hand side of eq. ( 15) respectively. For a piecewise linear polynomial basis, the number of degrees of freedom N DOF (associated to a finite element unknown in mesh Ω h ) is equal to the number of vertices N v of the mesh (N DOF = N v ). The solution vector x contains the finite element fields such that for a Ω h with N v vertices:

x =   w h 1,1 , w h 2,1 , w h 3,1 , w h 4,1 , . . . , w h 1,k , w h 2,k , w h 3,k , w h 4,k kth node , . . . , w h 1,Nv , w h 2,Nv , w h 3,Nv , w h 4,Nv    T ( 16
)
where {w h i,k } 4 i=1 are the values of the four finite element fields {w h i } 4 i=1 at the kth node of the mesh Ω h with

k = 1, 2, . . . , N v .
In the context of phase-field brittle fracture models, monolithic iterative methods (such as ours) may fail to perform in case of "brutal" crack evolution especially in the post-peak regimes [START_REF] Gerasimov | A line search assisted monolithic approach for phase-field computing of brittle fracture[END_REF]. Convergence can be attained by employing infinitesimally small continuation steps [START_REF] Alexander | The homotopy continuation method: numerically implementable topological procedures[END_REF] or by applying more advanced methods such as line-search assisted Newton iterations (Gerasimov and De Lorenzis, 2016), trust region methods [START_REF] Kopaničáková | A recursive multilevel trust region method with application to fully monolithic phase-field models of brittle fracture[END_REF], modified Newtons iterations [START_REF] Wick | Modified Newton methods for solving fully monolithic phase-field quasi-static brittle fracture propagation[END_REF], etc. For the chosen Picard iterations in this study we employ the infinitesimally small continuation steps in our numerical experiments.

Parallelization

For the vectorial finite element discretization procedure of the hybrid phase-field model, the linear system assembled from eq. ( 15) is of the form, Ax = b, with matrix A ∈ R m×m and vectors

x ∈ R m and b ∈ R m , with m = (n + 1) × N v .
Simulating real-world three-dimensional problems of fracture mechanics with phase-field models prerequisite fine meshes, consequently leading to a large number of finite element unknowns. As such, assembly and solving of Ax = b is associated with a large computational overhead in terms of time and memory usage.

To further complexify the situation, these linear systems need to assembled and solved multiple times (thousands) to account for the nonlinear and quasi-static iterations.

To subdue these problems associated with the computational overhead, we have implemented a message passing interface (MPI)-based domain decomposition method of parallelization built on top of our vectorial FEM discretization. This enables the use of distributed-memory framework of computing to assemble and solve Ax = b in parallel, and largely reduces the computational expenses. Briefly, domain decomposition applied in this study can be divided into three main steps:

(i) for a numerical simulation to be carried out on N p number of of processes, we start with partitioning of the global mesh Ω h into N p number of smaller meshes

{Ω h i } Np i=1 (subdomains) with more or less equal number of mesh nodes N v,i ≈ N v /N p . The subdomains {Ω h i } Np i=1
are then assigned to N p number of MPI processes (computational units). These {Ω h i } Np i=1 contain a single layer of ghost-nodes that is used to build MPI communication patterns between neighboring subdomains;

(ii) concurrently, within each MPI process, using eq. ( 15) with mesh Ω h i , we assemble a local matrix

A i ∈ R mi×mi and a local vector b i ∈ R mi , with m i = (n + 1) × N v,i .
As such, we create parallel tables of degrees of freedom with local and global views. Globally, N p number of matrices A i and vectors b i correspond to row-wise 1D partitioning [START_REF] Grama | Introduction to parallel computing[END_REF] of the system Ax = b. For a parallel simulation with large number of processes N p , since

{Ω h i } Np i=1
Ω h , the local matrices and vectors (A i ,b i ) are of smaller dimensions than the global (A,b). This implies N p times faster linear system assembly phase. This phase of assembly is 'embarrassingly-parallel', hence requires no communication between the involved MPI processes;

(iii) finally, by using N p MPI processes, the global problem is solved iteratively by the CG method which uses parallel matrix-vector products with a small MPI communication overhead.

Solution algorithm and preconditioning

The matrix A resulting from the vectorial FEM discretization introduced above is sparse, symmetric, and positive definite. As such, the Krylov subspace-based CG method -developed by [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] -makes a standard choice of iterative solution algorithm for solving the assembled linear system Ax = b. Although the CG method is well established for solving such kind of linear systems, it is likely to suffer from slow convergence for practical problems of interest with millions of degrees of freedom, which are expected for hybrid phase-field simulations of brittle fracture mechanics. To speed up the convergence of CG, preconditioning is used [START_REF] Málek | Preconditioning and the conjugate gradient method in the context of solving PDEs[END_REF].

Preconditioning here refers to transforming the assembled linear system Ax = b into an alternate linear system with favorable properties for the CG iterative solution. Generally speaking, preconditioning attempts to improve the spectral properties of A, and speeds up the convergence of CG by reducing the spectral condition number of the problem. Hence, preconditioning transforms the original linear system Ax = b into an alternate one with the same solution x, but which is easier to solve.

At the linear system level, preconditioning involves in solving:

M -1 Ax = M -1 b, (17)
where matrix M is the preconditioner associated to A. For successful preconditioning, the cond 2 (M -1 A) is expected to be less than the cond 2 (A) and/or the eigenvalues of M -1 A are clustered around 1. A perfect preconditioner would then be the case where M -1 = A -1 : in this case, the associated condition number cond 2 (M -1 A) is equal to one. However such preconditioner is not practical to use. Note that eq. ( 17), referring to the preconditioned linear system, represents a technique called the left-preconditioning. Besides this it is also possible to set up right-preconditioning:

AM -1 y = b, x = M -1 y, (18)
or to set up split-preconditioning:

M -1 1 AM -1 2 y = M -1 1 b, x = M -1 2 y, (19)
with preconditioner M = M 1 M 2 . Since the matrices M -1 A, AM -1 , and M -1 1 AM -1 2 have the same eigenvalues, therefore, the convergence of the CG method will be same for these different cases of preconditioners, eqs. ( 17) to ( 19) [START_REF] Benzi | Preconditioning techniques for large linear systems: a survey[END_REF]. As such, only left-preconditioning, eq. ( 17), is discussed and applied in this article. Also, for efficiently solving the linear systems of the hybrid phase-field model for fracture on parallel computing environment, the preconditioned CG method is a desirable solution technique since its main components -matrix-vector products, dot products, and norms -readily map to distributed memory implementations.

Broadly, preconditioning (constructing a preconditioner) can be achieved via two different approaches, the problem-specific approach and the purely algebraic approach. For a given problem, the problem-specific approach constructs the preconditioners that are based on approximated physics or method which is easier to solve but gives 'nearly' same solution than the given problem in hand, or preconditioners are constructed based on lower-order discretizations of the considered partial differential equations. On the other hand, the purely algebraic approach of preconditioning involves constructing the preconditioner solely by the information contained in coefficients of the assembled matrix A. These are often constructed after matrix A has been assembled. While it would be interesting to find/apply optimal preconditioners for fracture mechanics problems via the problem-specific approach, however such a preconditioner does not exist yet.

Although, one should note that the geometric multigrid preconditioner can be considered to be of the problem-specific kind and these have been used alongside the CG method for resolving fracture mechanics problems with phase-field models in the recent past (see e.g., [START_REF] Bilgen | A phase-field approach to conchoidal fracture[END_REF], [START_REF] Jodlbauer | Matrix-free multigrid solvers for phase-field fracture problems[END_REF]).

The problem-specific preconditioning is generally very sensitive to the details of the involved problem, and even modest changes in the problem can compromise the effectiveness of the solver. For such reasons, purely algebraic preconditioners are often preferred. Besides this, other advantages of the purely algebraic preconditioners over the problem-specific ones are discussed in [START_REF] Benzi | Preconditioning techniques for large linear systems: a survey[END_REF].

In this study, different preconditioners (purely algebraic ones provided in PETSc package Balay et al. ( 2019))

are studied and tested. In particular, attention is focused on:

(i) Jacobi preconditioning. Also known as the diagonal preconditioner, the Jacobi preconditioning is considered to be one of the simplest preconditioning techniques which involves setting M as the diagonal of A, such that M = diag(A 11 , . . . , A jj , . . . , A mm ), with A jj denoting the diagonal coefficients of A and j = 1, 2, . . . , m. For the parallel variant of Jacobi preconditioning applied in this study, the distributed M i ∈ R mi is applicable to A i so that M i = diag(A 11 , . . . , A jj , . . . , A mimi ), with i = 1, 2, . . . , N p , m = Np i=1 m i , and j = 1, 2, . . . , m i . This diagonal matrix M i is trivial to invert, as it simply involves replacing the coefficients of the matrix with its reciprocal,

M -1 i = {1/A jj } mi j=1 .
As such, a Jacobi preconditioner is computationally cheap and straightforward to construct. As discussed earlier, a reasonable preconditioner M is such that M -1 ≈ A -1 and is also cheap to construct. While the Jacobi preconditioning fulfills the latter condition, it is unlikely for a physics like the fracture mechanics that the M -1 constructed solely from the diagonal of A is a good approximate of A -1 .

(ii) Block Jacobi preconditioning. A computationally costlier but more robust preconditioner than the Jacobi preconditioner is the block Jacobi preconditioner. It involves setting M i ∈ R mi×mi to the diagonal block matrix of A i . The preconditioning matrix M i is defined by

M i = {A jk } mi,mi j=1,k=1 .
Such a preconditioner is particularly more effective than the Jacobi preconditioner since the block diagonals succeed in reflecting the nonzero structure of the coefficient matrix A hence hold considerable information of the underlying physics. As such, M -1 from block Jacobi preconditioning is a good approximate of A -1 in comparison to the Jacobi preconditioning. This statement should hold true for the linear systems assembled with the vectorial FEM discretization of the hybrid phase-field model for fracture. In this paper, to extract M -1 i from the local diagonal block matrices M i two options are used, namely the incomplete LU (ILU) factorization [START_REF] Dupont | An approximate factorization procedure for solving self-adjoint elliptic difference equations[END_REF]) and the incomplete Cholesky (IC) factorization [START_REF] Kershaw | The incomplete Cholesky-conjugate gradient method for the iterative solution of systems of linear equations[END_REF], both with zero fill-in. For the block Jacobi setup with ILU, the preconditioner M i = Li Ũi , where Li and Ũi are incomplete lower and upper triangular factors of local diagonal block matrix of A i , and M -1 i = ( Li Ũi ) -1 . And for the block Jacobi setup with IC, the preconditioner

M i = Li LT i , and M -1 i = ( Li LT i ) -1 .
(iii) Algebraic multigrid preconditioning. The above discussed preconditioning options are from the class of 'one-level' preconditioners; the AMG preconditioning discussed here falls under the category of 'multilevel' preconditioners. The AMG preconditioning is widely used, due to its optimal computational cost and algorithmic scalability. It accelerates the convergence of the CG method by building a hierarchy of coarser grids operators (solutions) from the fine grid operator. In other words, a series of independent coarse grid systems are used to accelerate the convergence of the overall fine grid system Ax = b. Note that, the AMG does not require access to the geometric grid, "grid" within the context of AMG preconditioning refers to set of indices of the unknown variables within the linear system Ax = b. The hierarchy of grids are obtained for a problem, starting from the finest grid level A and reducing the number of unknowns to get the coarser levels, until the coarsest level is small enough to be solved with exact LU factorization. Generally, the coarsest level contains a small number of unknowns in comparison to the finer levels, thus the cost of performing the exact LU factorization for the coarsest level is negligible in comparison to the overall solving time for the finest level linear system. Based on how coarsening is achieved, AMG is classified as the classical AMG [START_REF] Ruge | Algebraic multigrid[END_REF], the agglomeration AMG [START_REF] Jones | AMGe based on element agglomeration[END_REF], and the smoothed aggregation AMG [START_REF] Vaněk | Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems[END_REF]. Due to symmetric positive definite nature of the assembled matrix A, the smoothed aggregation AMG is employed in this study which was originally developed in [START_REF] Vaněk | Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems[END_REF] for elliptic systems. A key component of smoothed aggregation AMG is its smoother, which make the underlying error on the fine grids smooth so that it can be approximated accurately on the coarser grids. In other words, the smoother eliminates errors associated to large eigenvalues of the linear system, while the coarse grid correction eliminates the remaining errors associated to the small eigenvalues. There are many well known smoothers that are suitable for parallel AMG, c.f. Baker et al. (2011a,b), among these we employ the Richardson's and the Chebyshev smoothers in this paper.

Results and discussions

This section presents details, results, and discussions for the series of numerical experiments performed in this study. Our primary motive of performing different numerical experiments is to analyze, tune, and propose an optimal preconditioner for the hybrid phase-field fracture mechanics vectorial FEM solver capable of handling large-scale problems with millions of unknowns. We end this section by applying the tuned preconditioner to solve a complex three-dimensional brittle fracture problem with more than 64 million unknowns.

Hardware and software specifications

Different open-source frameworks (domain-specific language, mathematical libraries, etc.) were employed in this study to develop and implement the monolithic vectorial FEM solver and consequently perform the various sequential/parallel numerical experiments that are to appear in this section. Details of these open-source frameworks are provided in table 1.

The numerical experiments in this study were performed by utilizing computational resources that range from an ordinary "desktop PC" to a supercomputer "Inti" hosted at TGCC-CEA3 , Bruyères-le-Châtel,

France. The specifications of these two computational resources are provided in table 2. The "desktop PC" was used for executing the numerical experiments appearing in section 5.2, and for the numerical experiments in sections 5.3 and 5.4 Skylake nodes of "Inti" were used. For all the numerical experiments to follow, the CG algorithm convergence is assumed to be reached when the relative residual norm of the linear system is lower than 10 -10 or when the absolute residual norm is lower than 10 -10 . Additionally, the nonlinear Picard iteration at a given t are stopped for tolerances lower than 10 -8 .

Before proceeding further, in order to assist the analysis of results in the upcoming subsections, we introduce the following nomenclature: (i) to annotate different CPU timings, the total solver time t total is computed as the sum of three contributions:

t total = t assembly + t setup + t solve , where t assembly denotes the CPU time spent in the linear system assembly phase (i.e., to construct A and b), t setup is the CPU time spent in setting up the preconditioner operator M -1 , and t solve the CPU time spent in performing the CG iterations till convergence. We also denote the total solution phase CPU time t solution :

t solution = t setup + t solve ;
(ii) CG iteration number is denoted by k, the total number of CG iterations to convergence (within one iteration of the global iterative solving, i.e., from iteration jth to j + 1th in eq. ( 15)) is denoted by k c , and k c denotes the total number of CG iteration during one full simulation (i.e. for executing the full phase-field simulation);

(iii) at kth CG iteration, r k denotes the unpreconditioned residual, such that r k = ||Ax k -b|| 2 , with x k being the solution of the linear system;

(iv) the CG residual reduction rate R r k is defined by:

R r k = kc k=1 (r k-1 /r k ) k c -1 .

Two-dimensional numerical experiments: solver validation and preconditioner performance assessment

Commonly used numerical tests from literature (see e.g., [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], [START_REF] Liu | Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model[END_REF], [START_REF] Jeong | Phase field modeling of crack propagation under combined shear and tensile loading with hybrid formulation[END_REF], [START_REF] Hirshikesh | A FEniCS implementation of the phase field method for quasi-static brittle fracture[END_REF] to cite but a few), the two-dimensional (2D) single-edge notched tensile and shear fracture tests, are considered as benchmark problems in this subsection. From here-forth the tensile test is referenced as test 1 in the text.

Problem setting

The domain of interest is an initially cracked square plate (x, y) ∈ Ω = [0 cm, 1 cm] 2 (fig. 1a). With an initial parameter κ is set to 1 • 10 -6 and l 0 is assumed equal to 2h, where h is the characteristic size of the mesh

Ω h .
The unstructured Delaunay (triangular) meshes generated with Gmsh are used for solving the finite element problem of test 1. To establish mesh convergence, test 1 has been solved multiple times by varying the level of mesh refinements, details of these meshes are provided in table 3. The hierarchy of mesh refinements were generated by dividing each triangle in Ω h into four equal triangles. As such in table 3, we observe that with every refinement, the mesh size h halves and the number of triangles quadruple. The initial crack fields for the three mesh refinements (visualized using damage-field d) are presented in fig. 2.

Solver validation

Besides being used in the upcoming subsection to assess the performance of different preconditioners, we also used the test 1 to cross-validate and compare our vectorial FEM solver (sequential and parallel) against benchmark solutions of this test available in the literature. In fig. 3, the top surface reaction force F y versus applied displacements is plotted for finest mesh level 3 (detailed in table 3) and compared to reference hybrid phase-field solution from [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] and anisotropic phase-field solution from Miehe et al. (2010a). The force-displacement curve obtained from our computation at the finest mesh level 3 is in good agreement with the two reference solutions. This simulation was executed using 10 processes on the desktop PC mentioned in table 2. The parMETIS partitioned mesh with 10 subdomains is presented in fig. 1c. To further attest the correctness of the proposed vectorial FEM solver, we compare the errors in computing the top surface maximum reaction-force max(F y ) obtained from our solver against two different reference solutions provided in [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] and Miehe et al. (2010a). The last two columns of table 3 enumerate these errors at different mesh levels. At finest mesh level 3, these errors decrease down to 2% and 1%, respectively in comparison to the hybrid and the anisotropic phase-field solutions. This reduction in error while decreasing the mesh size suggests mesh convergence of vectorial FEM solutions towards the reference solutions from Ambati et al. ( 2014) and Miehe et al. (2010a).

Alongside the plot in fig. 3, four instantaneous snapshots of the calculated damage-fields are presented. These damage-fields are obtained from the simulation of test 1 at the finest mesh level 3. Damage-field evolution, crack initiation, and propagation can be observed in these snapshots. As expected, under extreme tensile loading, the crack can be seen to travel along a (almost) straight line dividing the square specimen into two (almost) equal halves.

Further validation is provided by performing the shear test for the same specimen of test 1. Now, the specimen is loaded in shear mode, i.e, the displacement Dirichlet condition (now on u 1 ) on the top edge is applied with an increment of ∆ū 1 = 1 • 10 -5 mm (see fig. 4a). Note that for this test and the previous one small displacement increments (< 1 • 10 -5 mm) are applied. Such small displacement increments are taken in view of easily achieving convergence for the nonlinear Picard iterations. For this test, the applied shear loading causes the top half of the specimen (above the initial crack) to be in compressive stress state while its bottom half (below the initial crack) remains in tensile stress state. As such, the crack is expected to propagate within the tensile stress state region (underneath the initial crack). Such a cracking trajectory obtained from our solver is compared to other reference cracking trajectories from literature in fig. 4b. The finest mesh level 3 was used to obtain this cracking trajectory, and the simulation was executed using 10 processes on the desktop PC mentioned in table 2.

In fig. 5, the force-displacements plot for the shear test at the finest mesh level 3 is presented and com- Note that for additional validation, other literature comparative tests (mode I, mode II, and mode III fracture) were also performed but these are not shown here for the sake of conciseness.

Effect of preconditioning on iteration count, condition number, and problem size

In the phase-field model, the coefficients of A are approximately dependent on (1 -d) 2 E (with E denoting the fourth-order elastic stiffness tensor, see eq. ( 1). As d approaches to 1 at certain (cracked) mesh nodes, this causes the condition number, cond 2 (A), to change during the numerical approximation of cracking processes. A direct effect of this change in the cond 2 (A) is sensed by the CG method, where the total number of iterations to converge (k c ) drastically change during the progression of the simulation.

Besides improving the convergence of the CG method, one of the roles of preconditioning operator M -1 is also to minimize the drastic change in k c during the full phase-field simulation, or in other words control the drastic change in cond 2 (A). In our numerical experiments from this subsection, we use such quality of the preconditioner as a metric to judge the efficiency of the preconditioning operator. Iterations. Figure 6 plots the complete evolution of k c observed while performing the full phase-field simulation, test 1 with level 1 mesh. The simulation was performed in parallel using 10 processes (MPI ranks) on the desktop PC from table 2. One can observe the effectiveness of preconditioning by noticing how k c gets reduced when CG method is used alongside different preconditioners.

In particular, k c with a simple Jacobi preconditioner, ranging from 593 to 920 reduces to a range of 152 to 202 when using a more robust block Jacobi preconditioner. Compared to the Jacobi and block Jacobi preconditioners, further drop in k c range can be observed with the AMG preconditioner, the k c iteration ranging now from 75 to 91.

For this test, the diagonal block matrices {M i } Np i=1 of the block Jacobi preconditioner were factorized by using zero fill-in incomplete Cholesky (IC) decomposition. For conciseness we shall drop the words "zero fill-in" and refer this preconditioner in the text as "block Jacobi IC". Further, the AMG preconditioner was in fact the smoothed aggregate AMG preconditioner [START_REF] Vaněk | Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems[END_REF] with the Chebyshev smoother and block Jacobi IC inner solver for the coarse grids and the coarsest grid being factorized by complete LU. The reason for choosing such kind of preconditioners will be discussed in upcoming section 5.3.

Within fig. 6, a strong rise in CG iterations to converge is observed as the load passes from 5.3 × 10 -3 to 6 × 10 -3 , i.e., during the damage propagation phase (illustrated in gray color in fig. 6). The strong rise in CG iteration count k c serves as an indicator that the conditioning of the assembled preconditioned linear systems (cond 2 (M -1 A)) must be changing drastically in the cracking zone.

In fig. 6, we also provide the data pertaining to the total CG iterations k c needed to perform the complete phase-field simulations (1490 pseudo-time steps) with Jacobi, block Jacobi, and AMG preconditioners.

Compared to the ∼ 1.1 million CG Jacobi iterations needed to completely execute the simulation, the block Jacobi (resp. AMG) completed the simulation in ∼ 1/4th (resp. ∼ 1/10th) number of Jacobi iterations.

Condition numbers. The speed of convergence for the CG method can be predicted by the knowledge of the condition number of the solved matrix A denoted by cond 2 (A) ≥ 1. The convergence of the CG method relies on the eigenvalue distribution of A, and cond 2 (A) provides a criterion for assessing the width of its spectrum. Following [START_REF] Basermann | Preconditioned CG methods for sparse matrices on massively parallel machines[END_REF], let us denote γ := ( cond 2 (A) -1)/( cond 2 (A) + 1). The distance between the CG solution x k at kth iterate from exact solution x is bounded:

x k -x 2 ≤ 2γ k cond 2 (A) x 0 -x 2 .
As such, the lower the condition number the faster the convergence of CG, or in other words lower k c . Now, to investigate on the effect of preconditioning operators, fig. 7 presents the condition number evolution during the full phase-field simulations of test 1 with mesh level 1. The cond 2 (M -1 A) is calculated by employing the Lanczos iteration algorithm for eigenvalue computation [START_REF] Balay | PETSc Users Manual[END_REF].

In practice, at each pseudo-time step, the Lanczos algorithm is used to calculate the eigenspectrum {λ i } Neigen i=1 of matrix M -1 A, with N eigen denoting the total number of eigenvalues calculated. Then, we use the

approximation: cond 2 (M -1 A) ≈ max|λ i | min|λ i |. Note that the eigenspectrum {λ i } Neigen i=1
provided by the Lanczos algorithm cannot be used for complete eigenvalue analysis; however it is intended here to assist in understanding the convergence behavior of the CG and to quantify the efficiency of the preconditioning operators [START_REF] Balay | PETSc Users Manual[END_REF][START_REF] Badri | Preconditioned krylov subspace methods for solving radiative transfer problems with scattering and reflection[END_REF].

Complementary to the results from fig. 6, the condition number evolution provided in fig. 7 reveals that AMG preconditioning is able to reduce the condition number deterioration during cracking. Hence the problem is easier to solve when compared with other preconditioning options. As suspected by the rise in k c from fig. 6, the plots from fig. 7 illustrates that cond 2 (M -1 A) rises exponentially as the crack starts to propagate (highlighted in gray color). Further analysis of these data reveals that in comparison to the Jacobi preconditioner, the average cond 2 (M -1 A) for the AMG (resp. the block Jacobi) is 1.8 (resp. 1.3) times smaller. This proves the superiority of the AMG preconditioner to "easily" solve the brittle fracture problems arising in the phase-field modeling. The main aim of this subsection was to study condition number degradation caused by cracking, while the condition numbers were the ones estimated from the monolithic linear system assembled for the hybrid phasefield model. It should be noted that the magnitude deference between displacement and damage fields partly contributes to poor conditioning of the studied linear systems. Dealing with a dimensionless formulation (i.e., scaling the problem) could lead to linear systems with relatively different condition numbers (possibly lower). In this study we do not deal with a dimensionless formulation, however for such cases as well cracking should lead to rise in condition number (relative to the non-cracked state). As such, the conclusion that the AMG preconditioner for CG is relatively superior to the Jacobi or the block Jacobi preconditioner should apply to linear systems assembled from a dimensionless formulations too.

Influence of the mesh refinement. In order to gauge the sensitivity of the three different preconditioning operators with respect to the refining of the mesh (problem size), table 4 presents total CG iterations to converge k c for a single pseudo-time step (9th) of the phase-field simulation test 1 with different mesh refinements. Since the highest order derivative for the equations of the phase-field model (eqs. ( 1) and ( 2)) looks like the Laplacian, as such, the condition number with the CG method will be O(h). So k c should increase as the square root of the problem size (number of elements).

Denoting the change in iteration count between simulations at mesh levels 1 and i by Dk c,1→i :

Dk c,1→i = k c,i /k c,1 = N e,i N e,1 ,
with k c,i and N e,i corresponding respectively to CG iterations to converge and number of elements for ith mesh level. In our numerical experiment, the mesh elements are quadrupled between mesh levels 1 and 2 (see table 3), therefore one observes approximately twice ( 65536/16384 = √ 4 = 2) increase in iterations Dk c,1→2 for the CG Jacobi method, since it is a weak preconditioner (see table 4 second row third column).

This trend continues with the other mesh levels for the Jacobi preconditioner. It indicates that the Jacobi preconditioner is a weak preconditioner for the phase-field problems of fracture in the sense of not being able to control the rise in CG iteration count Dk c,1→i between mesh refinements.

In contrast to what was observed for the Jacobi preconditioner, in table 4, we observe that the CG setup with the AMG preconditioning remains the least sensitive to mesh refinement. At the finest mesh refinement (level 5), for the AMG preconditioner we observe a Dk c,1→5 = 10.63 in comparison to Dk c,1→5 = 14.43 for the block Jacobi preconditioner, Dk c,1→5 = 15.65 for the Jacobi preconditioner, and Dk c,1→5 = 16

(theoretical) for no preconditioning. These results again highlight the superiority of the AMG preconditioner in comparison to the Jacobi or the block Jacobi preconditioners. Overall these results points out that for the AMG preconditioner CG convergence is least effected by mesh refinements. Although only a single pseudotime step (9th one) was used here to showcase the effect of mesh refinement on different preconditioners, from the results presented in figs. 6 and 7, we know that similar conclusions should apply for other pseudo-time steps of the full phase-field simulation.

Overall in this subsection, the preconditioning analysis -based on CG iterations to converge k c , total CG iterations needed for the full phase-field simulation k c , sensitivity of k c with mesh refinement, and condition number analysis -points towards AMG being an effective choice of preconditioner for solving brittle fracture problems using hybrid phase-field models in a vectorial FEM framework.

Three-dimensional numerical experiments: preconditioner tuning and scaling assessment

This subsection aims to test, tune, and compare the previously discussed preconditioning options for solving large-scale brittle fracture problems with the CG algorithm. In particular, two three-dimensional (3D)

problems with 32,764,068 (∼ 33 M) and 80,570,308 (∼ 81 M) unknowns are used as numerical tests in this subsection.

Problem setting. The considered problem in here is a 3D simulation of tensile damage in a pre-cracked specimen (fig. 8a) obtained by a unit extrusion (along z direction) of the 2D domain considered in test 1.

Similar to test 1, the bottom and top surfaces (xz-planes at y = 0 and y = 1) of the cubic geometry are respectively constrained and vertically loaded till failure. For the constrained plane ∂Ω D (x, y : y = 0, z), the displacements ū1 = ū2 = ū3 = 0; and for loading plane ∂Ω D (x, y : y = 1, z), the displacement Dirichlet condition is applied on u 2 with an increment ∆ū 2 = 1 • 10 -5 mm.

The unstructured finite element mesh used for the ∼ 33 M (resp. ∼ 81 M) contained 8,191,017 (resp. 20,142,577) nodes and 47,759,360 (resp. 118,308,864) tetrahedral elements. From here-forth in the text, we denote the ∼ 33 M and ∼ 81 M problems by test 2A and test 2B, respectively. Figure 8b shows the partitioned mesh used for the test 2A. This partitioned mesh contains 96 subdomains and is used for parallel numerical simulations with 96 processes.

Preconditioner tuning

To identify the optimal preconditioning option for large-scale phase-field fracture simulations, the linear systems obtained by the vectorial FEM discretization of tests 2A and 2B were solved by using the CG method setup with Jacobi, block Jacobi, and AMG preconditioners. Each simulations was carried out using 384 processes on 8 Skylake nodes of the Inti supercomputer. CPU timings. In fig. 9, comparative CPU timings for the different preconditioning options are presented in the form of stacked bar plots. Each stacked bar within the plot represents t total and has been color-coded with other CPU timings t assembly , t setup , and t solve . The color-coding of these stacked bars help in analyzing the CPU cost of different phases of the full finite element solution. Note that for performance/tuning/scaling analysis of preconditioners presented in this subsection we shall only be using a single pseudo-time step (8th one) from the full brittle fracture mechanics simulations. Further, to avoid any hardware related discrepancies, all timings are averaged based on three simulation runs for a specific case. Additionally, fig. 9 also provides information of total CG iterations to convergence k c for each preconditioning option.

Among the tested preconditioners, as expected, the CG setup with the Jacobi preconditioner yields high time to solution t total and high number of CG iterations to converge k c . Overall, the total solution time t total of 22.93 (resp. 65.41) seconds was observed for the numerical simulation of test 2A (resp. test 2B) solved parallely using 384 processes. This t total is decomposed into t assembly = 7.62 (resp. 18.98) seconds and t solution = 15.31 (resp. 46.43) seconds for the test 2A (resp. test 2B). Further, the Jacobi preconditioned CG took 561 (resp. 690) iterations to convergence for test 2A (resp. test 2B). Observe that the preconditioner setup cost t setup is negligible (∼ 1 % of the total solving time t total ). As the Jacobi preconditioning simply involves the use of diag(A) as the preconditioning operator M -1 , t setup is expected to be low. Figure 9: stacked bar plots for CPU timings highlighting the effect of different preconditioners on overall solver timings for two three-dimensional brittle fracture phase-field problems test 2A and test 2B solved using 384 processes. In the figure "BJacobi" stand for block Jacobi. The numbers in the brackets next to the bar plots represent overall solver timings ttotal in seconds (black) and the number of CG iterations kc to convergence (blue).

The CG setup with the block Jacobi preconditioner performed faster (less t solution and k c ) than the CG setup with the Jacobi preconditioner. Faster solution timings t solution were observed when the 384 inner diagonal block matrices {M i } 384 i=1 of the block Jacobi preconditioner were solved using the IC factorization, in comparison to the zero fill-in incomplete LU factorization (ILU) being used as the block diagonal matrix solver.

Further, the performance of the block Jacobi preconditioners can be improved by increasing the fill-in level ("n") of the ILU(n) or IC(n) used as the block diagonal solver to obtain M -1 . Note that with increase in n, ILU(n) (resp. IC(n)) approaches to complete LU (resp. Cholesky) factorization leading to a perfect preconditioner. However, since the goal of this article is to propose a preconditioner with good parallel performance for large-scale phase-field problems, we refrain from using the complete LU (resp. the complete Cholesky) factorization as diagonal block solvers to the block Jacobi preconditioner, and only use the zero fill-in versions of ILU and IC.

In comparison to the baseline run results from the Jacobi preconditioner, we observe that the t solution of block Jacobi ILU is 1.33 (resp. 1.25) times lower than the t solution of the Jacobi preconditioner for test 2A (resp. test 2B). We also observe that the t solution of block Jacobi IC is 1.61 (resp. 1.31) times lower than the Jacobi preconditioner for test 2A (resp. test 2B).

The CG setup with the AMG preconditioner, when properly tuned, proves to be the fastest preconditioner among the other preconditioning options tested in this article. This statement is asserted by values of total solving time t total and the CG iterations to k c for the AMG preconditioner provided in fig. 9. Moreover, as the combined displacement-damage linear systems (from vectorial FEM) to be solved are symmetric, the smoothed aggregate AMG preconditioner outruns (both in t total and k c ) the non-smoothed version (see table 5).

Among the smoothed aggregate, the Chebyshev smoother performed better (time-wise and iteration-wise) than the Richardson's smoother. To elaborate on this, in comparison to Richardson's smoother, the Chebyshev smoother uses the eigenvalue estimates for the smoothing procedure. Thanks to positive definite symmetric matrices provided by the vectorial FEM discretization, the smallest eigenvalues are easy to approximate by using CG as the eigenvalue estimator. As such, the Chebyshev smoother outperforms the Richardson's one.

Further, we tested which Krylov solver among the Jacobi, block Jacobi, and SOR performs the best as a coarse grid solver for the AMG preconditioner. Among them, the block Jacobi IC diagonal block solving outran the other two options. To further fine tune the AMG preconditioner, the coarse grid Krylov solver iterations were restricted to 1, 2, 3 , 4, and 5. It was revealed that 2 coarse grid solver iterations produce the fastest timings. In comparison to the baseline run results from the Jacobi preconditioner, we observe that the t solution of the tuned AMG preconditioner is 2.0 (resp. 1.9) times lower than the t solution of the Jacobi preconditioner for test 2B (resp. test 2A). Observe that in contrast to Jacobi preconditioner, where t setup was negligible in comparison to t solution , the same is not true for the AMG preconditioners. We observe almost 50 % of t solution consists of the preconditioner setup cost, t setup ≈ 0.5 t solution .

The results provided for the tests 2A and 2B in this subsection identify the smoothed aggregate AMG with coarse grid block Jacobi IC as the optimal preconditioning option for solving the monolithic (vectorial) finite element linear systems that arises from the considered hybrid phase-field model. In comparison to other preconditioners that were tested in this subsection, the tuned AMG preconditioner exhibited low iteration counts, low solution timings, and faster convergence rates. Note that from here on in the text, we shall refer this tuned AMG preconditioner by simply "AMG preconditioner".

Convergence history. In fig. 10, we show the convergence history for the three different preconditioners applied to CG as it approaches convergence at k c th iteration. Particularly, the relative residual r k is plotted as a function of the CG iteration number k ∈ [0, k c ] on a semilogarithmic plot. To obtain these plots, and for demonstration purposes, we simply used the test 2A solved on 384 processes. A nearly monotonic drop in residual can be observed for the AMG preconditioned CG solver, while the same is not true for the Jacobi or the block Jacobi preconditioned CG solvers. The residual reduction rate of the AMG is clearly superior to that of the Jacobi or the block Jacobi preconditioned CG solvers. Similar trends were observed (not shown here) for the test 2B.

For both problems, we also approximate the residual reduction rate R r k , which quantifies how rapidly a preconditioner converges. For the test 2A (resp. test 2B), R r k was observed to be equal to 2.09 (resp. 2.09) with the AMG preconditioner, whereas it was observed equal to be 1. 04 (resp. 1.03) and1 .11 (resp. 1.09) for the Jacobi and the block Jacobi preconditioners respectively. As such, with the AMG preconditioner, irrespective of the problem size, the residual approximately halves with each CG iterate as CG approaches to convergence. This is certainly a superior characteristic of the AMG preconditioner in caparison to the Jacobi and the block Jacobi preconditioners.

Scaling analysis

To analyze the parallelization efficiency of the proposed preconditioners, we now compare the strong parallel scaling performance of the Jacobi, block Jacobi IC, and tuned AMG preconditioners. The test 2A and 2B from the previous subsection, are used for the parallel scaling experiments which appear in this subsection.

However, to avoid confusion, these tests are renamed as tests 3A and 3B. To perform strong scaling analysis, 96 to 1,008 processes (MPI ranks) on Inti supercomputer (see table 2) were used. On this machine, test runs with 96 processes correspond to 2 Skylake nodes (N p = 2 × 48 = 96); similarly, 1,008 processes correspond to 21 Skylake nodes (N p = 21 × 48 = 1008).

In fig. 11, CPU timings data provided by the scaling analysis are presented in the form of stacked bar plots.

In particular, the CPU timings for linear system assembly and solving are highlighted. Additionally, the observed scaling for the different preconditioners is compared against the expected ideal scaling (theoretical).

Within the plots, the ideal parallel scaling efficiency and the observed one are also highlighted.

The scaling data suggests that all the considered preconditioning options lead to quasi-linear scaling behaviors when problem size is sufficiently large.

For test 3A, we observe excellent scaling efficiencies (∼ 95 %) for the three tested preconditioners, when the number of processes is between 96 to 384. However, a slight decrease in parallel efficiency was observed when processor count increased above 384. The following scaling efficiencies were observed for N p =1,008:

∼ 91 % for the CG Jacobi, ∼ 87 % for the CG block Jacobi IC, and ∼ 77 % for the CG AMG. In comparison to the parallel efficiency drop observed for the CG Jacobi, the efficiency drop observed for the CG AMG (resp. block Jacobi IC) preconditioner was significantly higher. This efficiency decrease is most likely caused by the higher MPI communication overhead of the CG AMG (resp. block Jacobi IC) solver in comparison to the CG Jacobi solver. Scaling-wise, all three preconditioning options have similar behavior. However, timing-wise, the CG AMG remained the fastest option. For instance, at N p =96, the simulation executed in 28.2, 37.0, and 57.8 seconds (solving time t solve ) with the CG AMG, the CG block Jacobi IC, and the CG Jacobi preconditioners, respectively. This makes the CG AMG solver 2.05 (resp. 1.31) times faster than the CG Jacobi (resp. block Jacobi IC) solver. At N p =1,008, the simulation solving times t solve of 4.1, 3.8, and 6.1 seconds were reported for the CG AMG, the CG block Jacobi IC, and the CG Jacobi solvers, respectively. In this case, due to higher overhead of the CG AMG for this particular problem, the CG block Jacobi IC is the better option among the three preconditioners.

As expected, in comparison to the test 3A, better strong scaling characteristics were observed for the test 3B for the three tested preconditioners. Note that, due to the large problem size, scaling data was measured in between 192 to 1,008 processes. The following scaling efficiencies were observed for N p =1,008: ∼ 99 % for the CG Jacobi, ∼ 97 % for the CG block Jacobi IC, and ∼ 96 % for the CG AMG. As observed for the test 3A, timing-wise, the CG AMG remained the fastest option among the three preconditioning options for the test 3B. For instance, at N p =196, the simulation solving times t solve of 47.4, 63.6, and 92.7 seconds were reported for the CG AMG, the CG block Jacobi IC, and the CG Jacobi, respectively. This makes the CG AMG 1.96 and 1.34 times faster than the CG Jacobi and the CG block Jacobi IC. In contrast to test 3A, at N p =1,008, the simulation solving times t solve of 9.2, 12.6, and 17.7 seconds were reported for the CG AMG, the CG block Jacobi IC, and the CG Jacobi, respectively, making the CG AMG 1.92 and 1.36 times faster than the CG Jacobi and the CG block Jacobi IC.

Overall, all three preconditioner tested in this subsection for CG showed excellent parallel strong scaling Figure 12: large-scale perforated medium test domain Ω and the partitioned mesh {Ω h i } 1008 i=1 . Within (a) the initial crack surface is highlighted in red color and the different sized holes can also be observed.

characteristics, however, timing-wise CG AMG remains the fastest option.

Large-scale application: a three-dimensional perforated medium

In this subsection, we present a mock application of perforated medium cracking. We demonstrate that by using the aforementioned parallel vectorial FEM formulation combined with the AMG preconditioner, we can "easily" perform the phase-field simulations of fracture propagation in a topologically complex medium, discretization of which leads to about 65,000,000 degrees of freedom.

Problem setting. A notched cubic specimen with the size of the side equal to 6 cm is submitted to tensile loading in the vertical direction (imposed via vertical displacements). The notch is placed at mid-height (y = 3 cm) of one side of the sample and is 1 cm in length. Its width is negligible if compared to the other dimensions of the computational domain.

To mimic a realistic perforated media, 44 spherical voids with different diameters are randomly arranged inside the specimen to induce a non-planar crack propagation. The generated geometry (fig. 12a) is such that the ratio between the solid phase and the voids is almost representative of concrete mortar. Mechanical parameters for performing the phase-field simulation are then chosen accordingly: λ = 1.94 kPa, µ = 2.45 kPa, and G c = 2.28 N mm -1 .

To induce damage initiation and evolution, the top surface of the numerical sample is loaded imposing a directly controlled vertical displacement increment of ∆ū 2 = 1 • 10 -3 mm, whereas the bottom face is fixed in all three directions ū1 = ū2 = ū3 = 0. Remaining lateral faces are stress-free. Similar to the preceding tests, these conditions serve as Dirichlet conditions for our finite element formulations.

The fine meshing constraints of this test case lead to a problem with 64,456,128 degrees of freedom. A finite element mesh with 16,114,032 nodes and 98, 295, 595 tetrahedral elements was used for solving this test with 1,008 processes on 21 Skylake nodes of Inti supercomputer. Figure 12b shows the partitioned mesh with 1,008 subdomains. To complete the full phase-field simulation, 865 total solves of the linear system were performed in less than 145 minutes (approximately two and a half hours). Considering that our solver possesses linear scaling characteristics, if this problem was solved on a sequential computer, hypothetically it would have taken (145 × 1008)/(60 × 24) ≈ 101 days to solve. This highlights the importance of parallel computations with effective preconditioning for solving large-scale realistic problems of fracture mechanics in reasonable times. Additionally, following the timing results from previous subsections, it is reasonable to say that CG setup with Jacobi and block Jacobi preconditioners would have taken more time and iterations to solve this particular test in comparison to the adopted AMG preconditioner.

Results. Without going into details of obtained results, briefly, the crack propagates from the notch to the right side of the specimen along the x axis. As a result, by the end of the simulation, the specimen is split into two halves. Figure 13 presents the displacement and the damage-fields after the crack is fully developed. Alongside these fields, the load-displacement plot of the top surface reaction force versus applied displacement has also been presented in fig. 13a. Typical tensile cracking behavior (mode I cracking) can be inferred from the curve, the reaction force peaks to a maximum value followed by rapid drop caused due to fast crack propagation. Notice from the displacement field (fig. 13b), for the whole upper half of the specimen u 2 ≈ 1 • 10 -3 (which is the imposed incremental displacement), while the lower half of the specimen remains fixed u 2 = 0. Complementary to this result, the damage-field shown in fig. 13c rightly suggests that the specimen is fully cracked. As such, we observe a free-body motion for the upper detached part of the specimen while the lower part remains static.

To highlight non-planar crack development, which was expected in this configuration, fig. 14 no more "planar", but deviates from the median plane due to the interaction with holes. Holes attract the crack due to local stress concentrations, thus leading to a quite realistic crack profile. Being able to trace such complex cracking trajectory infers robustness of our monolithic vectorial finite element procedure for solving the hybrid phase-field model.

Conclusions

In the current paper, we presented a monolithic finite element scheme of setting up the hybrid phase-field model for solving large-scale brittle fracture mechanics problems with millions of unknowns. To achieve a monolithic scheme, the combined discretization of the coupled damage-field equation and the augmented linear momentum equation was handled by employing the vectorial finite elements. All the finite element unknowns (displacement and damage fields) were consequently approximated via a single vectorial test and trial functions. Among others, one key advantage of using the vectorial finite element discretization for the phase-field model is that it is straightforward to parallelize via domain-decomposition method. As such in this study, we used distributed-memory computing for solving problems with extremely fine meshes which are a prerequisite for standard brittle fracture mechanics problems. Moreover, this study shows that when the vectorial finite element discretization is combined with the Picard iterations for handling nonlinearity posed by the coupled equations of the phase-field model, it leads to a single positive definite symmetric matrix system. This allows the use of Krylov subspace-based conjugate gradient method for solving (memory efficiently) these linear systems.

To improve convergence rates, consequently time to solution, of the conjugate gradient method applied to crack propagation problems with millions of unknowns, different preconditioning strategies were compared, analyzed, tuned, and discussed. Our analysis shows that the conjugate gradient method preconditioned with the tuned algebraic multigrid preconditioner is a robust and efficient scheme for solving the monolithic linear systems provided by the vectorial finite element discretization of the hybrid phase-field model for large-scale fracture mechanics problems.

By performing the condition number analysis on the monolithic linear systems for the brittle fracture mechanics problems, it was shown that crack initiation (during the phase-field simulations) causes the condition number to deteriorate (rise exponentially). This was consequently shown to impact the convergence behavior of the conjugate gradient solver, the iterations to converge grew drastically as the condition number to rose in the cracking zone. As such, the linear systems obtained for the hybrid phase-field problem within the cracking zone become more and more difficult to solve as the crack starts to move. Our numerical experiments prove that, in comparison to the conjugate gradient preconditioned with the Jacobi or block Jacobi preconditioners, the conjugate gradient preconditioned with the tuned algebraic multigrid preconditioner controls the rise in condition number during crack initiations. Hence, this latter preconditioner outruns the other two preconditioning options in the aspect of ease of solving a brittle fracture mechanics problem with the hybrid phase-field method. Further, numerical experiments also revealed that this tuned algebraic multigrid preconditioner remains least sensitive to mesh refinements. Considering that brittle fracture problems have a prerequisite of extremely fine meshing for replicating the cracking physics accurately, this quality makes the tune algebraic multigrid preconditioner an excellent choice over the Jacobi or block Jacobi preconditioners for the conjugate gradient method.

Results from two large-scale hybrid phase-field numerical simulations (with ∼ 32 and ∼ 81 million unknowns) performed for three-dimensional brittle fracture mechanics reveal that timing-wise, the tuned conjugate gradient algebraic multigrid preconditioner is approximately two (resp. 1.3) times faster than the conjugate gradient method setup with the Jacobi (resp. block Jacobi) preconditioner. Further, these large-scale tests were also used to show that the tuned algebraic multigrid preconditioner has the highest residual reduction rate (quality of rapidly decreasing the residual ||Ax -b|| to reach faster convergence). It was equally revealed that the algebraic multigrid preconditioner alongside the block Jacobi and Jacobi preconditioners possess excellent parallel strong scaling characteristics (more than 90 % efficiency) when these are applied to solve the large-scale phase-field fracture problems which involve linear systems with millions of degrees of freedom. However, timing-wise, the tuned algebraic multigrid preconditioner always remains the fastest option.

As an application case, a three-dimensional numerical simulation mimicking perforated media cracking which involved more than sixty-four million unknowns was "easily" solved by employing the vectorial finite element scheme. Notably, these are scale which would have otherwise not been possible to achieve if the tune preconditioning options for the conjugate gradient method was not used, which efficiently uses the high-performance computing resources.

  crack and a constrained bottom edge ∂Ω D (x, y : y = 0), the plate is subject to increasing displacements on its top edge ∂Ω D (x, y : y = 1) until the plate fully cracks open. The initial crack is placed at the center of the plate, i.e., ∂Ω D (x : 0 ≤ x ≤ 0.5, y : y = 0.5). These boundary conditions are also illustrated in fig.1a. The plate material is characterized by λ = 121.15 kPa, µ = 80.77 kPa, and G c = 2.7 kN mm -1 . Concerning the computational specifications of test 1, the displacement discontinuity imposed by the initial crack was modeled by nearly overlapping (tolerance δy = 10 -7 m) Dirichlet nodes placed along the cracks edge ∂Ω h D (x : 0 ≤ x ≤ 0.5, y : y = 0.5 ± δy) within Ω h . For illustration proposes, a coarse grid featuring Dirichlet nodes for the initial crack of test 1 is presented in fig. 1b. The displacement Dirichlet condition on the top edge is applied with an increment of ∆ū 2 = 1 • 10 -5 mm up to u 2 = 5 • 10 -3 mm and ∆ū 2 = 1 • 10 -6 mm up to failure of the specimen. For the lower edge, the constrained displacement Dirichlet conditions ū1 = ū2 = 0 are applied. Further, for test 1 and for all the simulations that appear in this study, u1 = u2 = 0 ∂Ω D : u2 = u2 + ∆u2 ∂Ω D : d = 1 (a) domain Ω. (b) a coarse mesh Ω h . (c) partitioned mesh {Ω h i } 10 i=1 .

Figure 1 :

 1 Figure 1: domain Ω, mesh Ω h , and partitioned mesh {Ω h i } 10 i=1 for test 1. (a) also illustrates the boundary conditions applied to test 1. (b) represents a coarse unstructured finite element mesh with 'nearly' duplicate Dirichlet nodes for the initial crack.

Figure 2 :

 2 Figure 2: initial crack visualization of test 1 via damage-field d at different mesh levels.

Figure 3

 3 Figure 3: force-displacement plot for the two-dimensional single-edge notched tensile fracture simulation, test 1. The orange solid line refers to our vectorial FEM solution, the square markers denote the reference hybrid phase-field solution from Ambati et al. (2014) and the triangular markers denote anisotropic phase-field solution from Miehe et al. (2010b).

Figure 4 :

 4 Figure 4: shear cracking test domain Ω and final cracking trajectory in compared to the trajectories obtained by the hybrid phase-field model[START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], the anisotropic phase-field model(Miehe et al., 2010a), and adaptive phase-field modeling[START_REF] Hirshikesh | Adaptive phase-field modeling of brittle fracture using the scaled boundary finite element method[END_REF].

  fig. 4b.

Figure 5

 5 Figure 5: force-displacement plot for the two-dimensional single-edge notched shear fracture simulation, test 1a. The orange solid line refers to our vectorial FEM solution, the square markers denote the reference hybrid phase-field solution from Ambati et al. (2014) and the triangular markers denote anisotropic phase-field solution from Miehe et al. (2010b).

Figure 6 :

 6 Figure6: evolution of the CG iterations to converge kc during execution of the full phase-field simulation (test 1 with level 1 mesh).

0 1 •

 1 Figure 7: evolution of condition number cond 2 (M -1 A) during execution of the full phase-field simulation (test 1 with level 1 mesh).

  partitioned mesh {Ω h i } 96 i=1 of test 2A.

Figure 8 :

 8 Figure 8: three-dimensional brittle cracking problem geometry and the partitioned finite element mesh {Ω h i } 96 i=1 . The precracked surfaces in (a) are highlighted in red color .

Figure 10 :

 10 Figure10: conjugate gradient method convergence history for the Jacobi, block Jacobi and algebraic multigrid preconditioners.

Figure 11 :

 11 Figure11: bar plots depicting the strong scaling characteristics of different preconditioners up-till 1,008 cores for the threedimensional brittle fracture problems, tests 3A and 3B, solved on Skylake nodes of Inti supercomputer. The numbers in the brackets next to the bar plots represent ideal speedup (black) and observed speedup (red). In the figure "BJacobi" stands for block Jacobi preconditioner.

  partitioned mesh {Ω h i } 1008 i=1 .

Figure 13 :

 13 Figure 13: results of the three-dimensional cracking of the perforated medium. Within the truncated visualizations (right ones in (b) and (c)) the spherical hole wire-frame meshes are colored according to the vertical displacement field u 2 .

  depicts the damage-field and crack surface development for three phases of the simulation, d is filtered for values close to unity in order to clearly identify the represented pseudo-crack. The color map in fig. 14-left helps to analyze the vertical non-planar deviation of this pseudo-crack (fracture deviation y d-d0 ). On the derived cracking surface, at a particular mesh node (P = P (x, y, z)), the fracture deviation y d-d0 ∈ R is simply obtained by y d-d0 = y -30 mm. Results from fig. 14-left suggests that due to the presence of holes, the crack deviates (compared to a planar crack) by approximately half a centimeter at certain places. As clearly illustrated by traversal cross-sections of fig. 14-right, despite what is obtained in the previous case, the damage band is

  pseudo time-step 256 (ū2 = 2.56 × 10 -1 mm). (b) pseudo time-step 261 (ū2 = 2.61 × 10 -1 mm).

  (c) pseudo time-step 266 (ū2 = 2.66 × 10 -1 mm).

Figure 14 :

 14 Figure 14: results depicting crack propagation within the three-dimensional perforated medium test. Left: truncated damagefield with fracture deviation color map; Right: different damage-field cross-sections (xy-planes) within the bulk of the material. For each (a), (b), and (c) the three xy-planes -top, middle, and bottom -are defined by top: (x, y : 20 > y > 40, z : z = 10), middle: (x, y : 20 > y > 40, z : z = 30), and bottom: (x, y : 20 > y > 40, z : z = 50).

Table 1 :

 1 packages used for setting up the monolithic vectorial FEM fracture mechanics solver. Note that, except for Gmsh all other packages were compiled with the Message Passing Interface standard (MPI) library[START_REF] Gropp | Using MPI: portable parallel programming with the message-passing interface[END_REF] to set up workload in distributed parallel computing framework.

	package	version operation

Table 2 :

 2 details of computational resources used for performing the serial/parallel numerical experiments in this study. In the table header "PU" denotes processing units (cores/threads).

	name	#nodes #PU architecture	test runs	interconnection
	desktop PC Inti	1 21	10 48	Intel Xeon E5-2680 v2-@2.80GHz sequential/parallel -Intel Xeon E5-8890 v4-@3.4GHz parallel InfiniBand EDR

Table 3 :

 3 characteristics and computational details for the different finite element meshes used for test 1. E1 max(Fy ) and E2 max(Fy ) are are the the maximum reaction force errors computed against references[START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] andMiehe et al.

	(2010a), respectively.
	mesh	nodes	triangles h	N DOF	Nnz	E1 max(f )	E2 max(f )
	level 1 8,353 level 2 33,089 level 3 131,713 262,144 16,384 65,536	0.0156 25,059 0.0078 99,267 0.0039 395,139 8,274,825 2.00% 520,425 15.62% 2,073,033 6.71%	14.78% 5.89% 1.01%

Table 4 :

 4 total CG iterations to converge for different preconditioning options. Numbers in brackets indicate the rise in iteration count Dk c,1→i . k Jacobi refer to the total number of CG iterations to converge kc for Jacobi, block Jacobi, and AMG preconditioners respectively. Subscript "BJacobi" denotes block Jacobi.

	c	, k B-Jacobi c	, and k AMG c
		level (i) N e,i /N e,1	k Jacobi c	k BJacobi c	k AMG c
		1 2 3 4 5	1 4 16 64 256	663 (1) 1,312 (1.98) 2,612 (3.94) 5,191 (7.83) 10,375 (15.65) 2,395 (14.43) 797 (10.63) 166 (1) 75 (1) 312 (1.88) 126 (1.68) 604 (3.64) 258 (3.45) 1,201 (7.24) 475 (6.34)

Table 5 :

 5 solution phase timings tsolution and CG iterations kc observed by varying different tuning parameters of the AMG preconditioner.

		max coarse space Krylov iterations	smoothed	smoother
	problem	1	2	3	4	5	yes	no	Chebyshev	Richardson
	Timings
	test 2A test 2B	20.2 62.2	7.8 24.5	15.4 77.1	18.6 83.9	20.4 95.7	7.8 24.5	35.5 108.4	7.8 24.5	31.0 97.3
	Iterations
	test 2A test 2B	176 214	42 48	55 108	34 89	30 84	42 48	155 187	42 48	112 214

We recall that, given a second order tensor A, its positive (A + ) and negative (A -) parts are defined as A ± = n i=1 A i ± e i ⊗ e i , where {A i } n i=1 and {e i } n i=1 denote the eigenvalues and eigenvectors of A, and symbol ⊗ indicates the dyadic product between vectors.

TGCC-CEA, "Très Grand Centre de Calcul" is an infrastructure for scientific high-performance computing and Big Data, able to host petascale supercomputers.

Acknowledgments

The authors acknowledge the financial support of the Cross-Disciplinary Program on Numerical Simulation of CEA, the French Alternative Energies and Atomic Energy Commission. We thank TGCC-CEA, Bruyères-Le-Châtel, France, for providing us compute time on the Inti supercomputer. We also thank Dr. Pierre Jolivet from IRIT, CNRS, Toulouse, France, for fruitful discussions on preconditioning/parallelization. The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.