Block preconditioning strategies for generalized continuum models with micropolar and nonlocal damage formulations

In this work, preconditioning strategies are developed in the context of generalized continuum formulations used to regularize multifield models for simulating localized failure of quasi‐brittle materials. Specifically, a micropolar continuum extended by a nonlocal damage formulation is considered for regularizing both, shear dominated failure and tensile cracking. For such models, additional microrotation and nonlocal damage fields, and their interactions, increase the complexity and size of the arising linear systems. This increases the demand for specialized preconditioning strategies when iterative solvers are adopted. Herein, a block preconditioning strategy, employing algebraic multigrid methods (AMG) for approximating the application of sub‐block inverses, is developed and tested in three steps. First, a block preconditioner is introduced for linear systems resulting from micropolar models. For this case, a simple sparse Schur complement approximation, which is practical to compute, is proposed and analyzed. It is tested for three different discretizations. Second, the developed preconditioner is extended to reflect the additional nonlocal damage field. This extended three‐field preconditioner is tested on the simulation of a compression test on a sandstone sample. All numerical tests show an improved performance of the block preconditioning approach in comparison to a black‐box monolithic AMG approach. Finally, a problem‐adapted preconditioner setup strategy is proposed, which involves a reuse of the multigrid hierarchy during nonlinear iterations, and additionally accounts for the different stages occurring in the simulation of localized failure. The problem‐adapted preconditioning strategy has the potential to further reduce the total computation time.


INTRODUCTION
Generalized continuum formulations, for which efficient solution techniques are studied in the present paper, are wellestablished approaches for the simulation of failure of cohesive frictional materials.This class of materials includes different types of solid geomaterials like rock and concrete, tough ceramics, but also granular materials like sands or, more general, soils.Depending on the specific material and the type of loading, failure can occur as shear failure, or in the form of quasi-brittle cracking.Generalized continuum models fall in the category of multifield approaches.They involve several coupled primary unknowns, that is, displacements and microrotations for micropolar formulations, 1,2 displacements and nonlocal damage variables for gradient-enhanced nonlocal models, 3,4 or displacements, nonlocal damage variables and microrotations in the context of combined micropolar/nonlocal damage approaches. 5,6Micromorphic continua, 7 as well as phase-field models 8 are also representatives of this category of multifield approaches.The focus of the present publication is, however, on the aforementioned micropolar and combined micropolar and gradient-enhanced nonlocal approaches.
For the considered generalized continuum models, the unknown fields are determined by the solution of a set of balance equations, among which are the balances of linear and angular momentum, and possibly a third balance equation of the continuum enhanced by gradients of internal variables for gradient-enhanced nonlocal models.The system of balance equations is closed by several kinematic and material-specific constitutive equations, resulting in a coupled set of partial differential equations.Due to the presence of large inelastic deformations, material and geometric nonlinearities are prominent in practical applications of the method.
Discretization of the corresponding coupled boundary value problem using the finite element method yields a set of nonlinear equations, which is commonly treated using a load incrementation procedure combined with Newton-Raphson solves in each load step.In this context, the most time-consuming step is solving the arising series of large and sparse linear systems, one in each step of each Newton iteration performed in the course of the simulation.Since the considered models involve several coupled partial differential equations and independent fields, the resulting linear equation systems are commonly characterized by a large number of degrees of freedom (DOF) and a distinct block structure.In this context, each diagonal sub-block in the system matrix is associated with a field equation of an uncoupled single field problem, while coupling between the fields is reflected in the off-diagonal coupling blocks.
The tremendous advances in the formulation, development, and application of generalized continuum models have driven the demand for problem-specific efficient and robust linear solvers.Due to memory requirements, iterative solvers are often preferred for large linear systems.Since multifield problems in general, and generalized continuum models in particular, are notorious for resulting in large linear systems, this is especially true for the models considered in the present publication.
The success of an iterative solution approach, however, is highly dependent on the quality of the preconditioner employed in the solution process.For instance, if a linear system is derived from a discrete Poisson-like equation, it is well known that the condition number depends on the discretization.This dependence directly affects the convergence properties of Krylov iterative solvers. 9,10The statement that preconditioning is crucial for successful iterative solutions of linear systems particularly applies to multifield problems.2][13] In computational mechanics, the preconditioning challenge appears for instance in contact problems where a saddle point structure emerges when Lagrange multipliers are used to model the contact constraints. 146][17] For these problems, the monolithic solution, that is, solving for all fields simultaneously, is often considered to avoid splitting errors. 17or the modeling of brittle fracture using phase-field approaches, research on preconditioning strategies is reported in a number of places in the literature.For instance, in Ref.,18 a block preconditioner is formulated for phase-field models.The preconditioner is based on an approximation of the Schur complement, and it is used within an outer MINRES Krylov solver.In Ref.,19 the performance of algebraic multigrid preconditioning for phase-field models is investigated.Work on preconditioning for three field problems arising in the context of phase-field modeling of fracture in incompressible materials can be found in Ref. 20 In Ref.,21 block-diagonal preconditioning is used in matrix-free approaches for the solution of phase-field fracture problems, and nonlinear preconditioning for phase-field modeling of fracture is proposed in Ref. 22  Tailored preconditioners for specific multifield problems therefore represent a vast area of research, and the present paper aims at contributing to this field of study for the special case of generalized continuum models.Since, as it will be shown by the numerical experiments included in the present paper, single-field preconditioners cannot be effectively applied to the multifield systems arising from generalized continuum models with micropolar and nonlocal damage formulations, different strategies have to be considered.However, before specific strategies are introduced and investigated in this study, basic notions for preconditioning of linear systems will be introduced in Section 1.1, important aspects of algebraic multigrid methods will be reviewed in Section 1.2, and the main governing balance equations for the considered generalized continuum models will be summarized in Section 1.3.Based on these considerations, an outline of the paper will be given in Section 1.4.

Preconditioning of linear systems
A well-chosen preconditioner is an essential step in ensuring the success of iterative solvers and enabling the practical computation of large-scale problems.The preconditioner transforms the original system into an alternative with favorable properties for the iterative solution.For left preconditioning, instead of computing the linear system involving a square  ×  matrix  ∈ ℝ × , the equation is solved, where  −1 ∈ ℝ × is the nonsingular preconditioner matrix,  ∈ ℝ  is a vector with the Newton corrections to the field variables, and  ∈ ℝ  is an out-of-balance vector.Ideally, the preconditioner  −1 should be an approximation for the inverse of the system matrix .
In practice, in the context of a Krylov iterative solver (e.g., GMRES 23 ), the preconditioner matrix  −1 is not explicitly formed.Instead, only its application to a vector is needed.For a known vector  ∈ ℝ  , the preconditioner action is defined by This is equivalent to solving a linear system,  = , for an unknown vector  ∈ ℝ  .This action, then, can be approximated by another iterative solver technique.One popular approach is to use algebraic multigrid methods, which are often optimal for elliptic problems. 24

General concepts of algebraic multigrid methods
Algebraic multigrid method (AMG) 25 is a popular iterative method for sparse linear systems.Originally, it was developed to solve discretized scalar elliptic differential equations.It can be used as a standalone solver or as a preconditioner for a Krylov subspace method.
In the present paper, the respective single-field preconditioner will be used as a reference method for evaluating the performance of the proposed preconditioning approaches.AMG is also a powerful method to compute approximate subblock inverses in block preconditioners.It, thus, forms the basis for the preconditioning approaches that will be developed within the present publication.
The main components of AMG are a relaxation step and a coarse grid correction.The relaxation step, also known as level smoother, can be, for instance, a few sweeps of classic fixed-point iterations within the framework of Jacobi or Gauss-Seidel methods.This step reduces specific error modes efficiently.The remaining error components, on the other hand, reduce slowly.Those are defined as algebraically smooth errors. 26The coarse grid correction in a complementary way reduces the smooth errors.

Overview
The method involves two phases, the setup phase of the multilevel 1 hierarchy and the solution phase.During the setup phase, the coarse level can be constructed based on a coarse-fine (CF) split of the system's DOF, 25 or by aggregation. 27In this work, we favor the aggregation strategy, which is part of the smoothed aggregation algorithm (SA-AMG) 2 .Considering a linear system  =  with ,  ∈ ℝ   , the setup output is the system matrix on all levels { () }, with  denoting the level and  (0) =  is the original matrix.On the subsequent levels, the system is progressively smaller.Also, as a setup output, the corresponding transfer operators for prolongation and restriction are generated for all coarser levels.After constructing the hierarchy, the solution process can be initiated.The process consists of a presmoothing on the fine level system with  1 relaxation iterations on the system  (0)  = .It results in an approximation  ( 1 ) of the solution.Then, the residual  (0) =  −  (0)  (𝜈 1 ) is transferred to the first coarse level by  (1) =  (1) with  () ∈ ℝ  +1 ×  as the restriction operator for level  and   the number of DOF at level .Notice that the restriction takes a vector from the fine level with dimension   to a coarse level with a smaller size  +1 .On this coarse level, the error-residual relation is used to form a coarse-level problem  (1) (1) =  (1) , where the error  (1) is computed or the new residual is transferred to another coarse level by recursion.Finally, when a coarse level result is obtained, it is interpolated to the fine level to correct the solution, which is then subjected to a postsmoother step.In a two-level scheme, the correction to the solution is given by with  ( 1 ,) representing the updated solution by the coarsening procedure and  () ∈ ℝ   × +1 the prolongator operator that interpolates from a coarse level to a fine level.Now, the prolongator takes a vector from the coarse level with dimension  +1 to a vector with larger dimension   , opposite of the restriction operator.Postsmoothing is performed by  2 iterations of a relaxation on the system  (0)  = , with initial guess  ( 1 ,) , which results in the final approximation after the multigrid cycle  ( 1 ,, 2 ) .

Setup phase
After this overview, we review some aspects specific to the setup phase, namely, the strength of connection and construction of all required transfer operators.In the coarsening process, a measure of connection strength defines which DOF should be grouped.A usual measure is based on the matrix entries.The next step of the setup after forming the aggregates is to construct the prolongator operator.This operator interpolates the error from the coarse to the fine level.The prolongator construction in the SA-AMG directly considers information about slowly converging errors.This information is related to the low-energy modes of the system, which in elasticity case is represented by the rigid body motions.Details about the significance of low-energy modes in AMG can be found in Refs.28, 29.After the prolongator operator is constructed, a few smoother iterations are applied to its columns, which is known to improve the convergence of the method. 27he restriction operator  can be constructed with the transpose of the prolongator operator  =   .This is known as the Galerkin approach, and it works well for symmetric positive definite problems.Although it can also be applied 1 In the algebraic multigrid method, the coarse grid has no relation with the original discretization grid, or mesh, hence the level denomination is more general and is used instead of grid. 2 The SA-AMG allows a natural inclusion of near null space information into the transfer operators.The near null space information is related to slowly converging error modes.This ensures that the coarse level correction is interpolated to the fine level directly on the smooth errors.See Ref. 27 for more details.
to nonsymmetric problems, it does not guarantee convergence. 26Alternatively, for nonsymmetric problems, a Petrov-Galerkin coarsening is adequate. 29In this case, the transfer operators are not equal to  ≠   .In our experiments, we opted for a Galerkin approach, which provided a good balance between performance and setup cost.

Generalized continuum models for material failure
As an example in this work, we consider the micropolar continuum, a generalization of the classical continuum considering an independent field of microrotations.Neglecting body forces and body couples for the sake of simplicity, it is expressed by a set of two coupled partial differential equations, that is, (i) the classical balance of linear momentum (ii) the balance of momentum moments with   denoting the macroscopic Cauchy stress tensor,   the couple stress tensor, and   denoting the Levi-Civita permutation symbol.
The independent fields are displacement and microrotation vectors,   and   , which are related to the stress tensors via proper constitutive models: =   ( , ,   ,  , ) .
Compared to the classical continuum, the involved stress tensors are nonsymmetric in general.
Due to the consideration of independent material rotations, micropolar continuum models are particularly well suited for modeling shear dominated failure, for example, Ref.,30 and hence received a certain amount of attention in the geomechanics community.However, micropolar continuum models are less suited for modeling fracture in tension or similar fracture modes, which do not involve rotations of the microstructure.
This deficiency can be resolved, for example, by considering an additional nonlocal scalar field αd describing the integrity of the material.Commonly, this scalar field is governed by an additional driven by a local damage parameter  d on the right-hand side, and with a material length scale parameter  d defining the radius of nonlocal interaction.
Then, this additional independent field, governed by a second-order PDE, is commonly used for gradient-enhanced (or mathematically similar phase field) fracture approaches, 31 considering material damage as a nonlocal process by introducing a dependence of the stress tensors on αd : =   ( , ,   ,  , , αd ) =   ( , ,   ,  , , αd ) .
The evolution of  d is governed by the local constitutive response of a material, discussed in later sections.
The combination of the gradient-enhanced continuum and the micropolar continuum was presented recently in Ref.32 in the context of a unified finite strain damage-plasticity framework.It was shown that this framework allows us to realistically represent material failure under various loading conditions in numerical simulations, including cracking in tension and shear dominated failure.
Clearly, solving the resulting equation system for such a unified approach in numerical simulations is challenging: In 3D simulations, each material particle is endowed with 7 DOF, and the resulting system is characterized by a 3 × 3 block structure.Even if, in a micropolar continuum, damage is not considered, each material particle is still endowed with 6 DOF, and the resulting system has a 2 × 2 block structure.This confirms the demand for improved precondition strategies tailored to this class of problems.

Paper outline
For solving the coupled equation systems resulting from generalized continuum models efficiently by means of iterative solvers, advanced preconditioning strategies are indispensable.In particular for a micropolar continuum with an independent field of microrotations, such strategies have not yet been explored systematically.This is the focus of the present work.Therefore, in a first step, preconditioning will be investigated in Section 2 for an elasto-plastic micropolar continuum featuring displacements and micro-rotations.The required constitutive equations will be summarized briefly, and a suitable block preconditioning strategy will be developed.Then, the performance of the proposed preconditioner will be evaluated in comparison to single-field preconditioning for the analysis of shear band formation in plane strain compression tests on samples of Gosford sandstone.
Then, in a second step, the focus turns to the unified micropolar gradient-enhanced continuum in Section 3.After a brief review of micropolar damage-plasticity modeling for sandstone with three independent fields, an extended block preconditioner for the unified micropolar gradient-enhanced continuum will be proposed and assessed, again for the chosen test case of a plane strain compression test on Gosford sandstone.
In addition, in Section 4, the developed approach for micropolar damage-plasticity modeling is refined further by proposing a new problem-adapted strategy for the preconditioner setup in the context of numerical simulations of material failure.The effectiveness of the method is demonstrated by further numerical results.Section 5 concludes with a summary and recommendations for future research activities.

DEVELOPMENT OF A PRECONDITIONING STRATEGY FOR MICROPOLAR CONTINUUM MODELS
This section describes the development of a block preconditioner tailored to the micropolar continuum model.Although the preconditioner will be tested also in a materially nonlinear, large deformation setting later on, it is developed based on considerations for the special case of small deformations only.Therefore, to begin with, basic constitutive and kinematic relations for the micropolar continuum are summarized for this restricted setting.The respective variational form of the governing equations in the geometrically linear setting will be briefly reviewed, including the results of its spatial discretization using the finite element method.A proper preconditioner is then deduced based on properties of the set of algebraic equations obtained in the linear elastic regime.The performance of the proposed preconditioner is then evaluated for plane strain compression of a sandstone sample, which is modeled in this section using an elastic-plastic micropolar approach without damage.The interested reader is referred to the paper 32 for a full description of the governing equations formulated in a geometrically exact setting.Subsequently, bold capital letters are used to designate some specific finite element constructs such as interpolation operators and vectors with nodal values.

Summary of a geometrically linear micropolar elasto-plasticity formulation
In the linear elastic regime, the two constitutive relations, ( 9) and (10), are expressed in terms of six material constants as with  and  as the Lamé parameters; , ,and  are the micropolar parameters; and   is the coupling modulus, which controls the level of asymmetry of the shear components of the stress tensor.The micropolar parameters can be arranged as derived constants, such as the characteristic length for bending Furthermore, by restricting the range of those parameters to exclude nonphysical material response, one ends up with the characteristic length for bending  b and the coupling modulus   as independent parameters for the micropolar elasticity.
The micropolar deformation is defined by the asymmetric strain tensor   and a twist, or curvature, tensor   .Under small displacement assumptions, those are defined as Plastic deformation occurs if a stress state reaches the zero-isosurface of a Drucker-Prager yield function   , which is formulated in terms of a generalized mean stress invariant   , a generalized deviatoric stress invariant , and using an equivalent plastic deformation measure   .
The definition of the generalized invariants and the plastic deformation measure can be found in Ref.5 for the geometrically linear setting, and in Ref. 32 for the more general large deformation setting.In both cases, the definition of the generalized deviatoric stress invariant  contains a further length parameter  J2 .The material parameters defining the shape of the yield surface are cohesion , friction angle , and the hardening variables ℎ Δ and ℎ  .Plastic flow is governed by a flow rule formulated with a nonassociated plastic potential   given as It depends on the dilation angle , and a material parameter ℎ  defining the change in the plastic potential with increasing plastic deformation.Again, the reader is referred to Refs.5, 32 for more detailed information on the employed flow rules in the geometrically linear and nonlinear setting.

The discrete system in the geometrically linear setting
The discrete system of equations is obtained with the finite element method.First, the balance Equations ( 7) and ( 8) are recast in the weak form with   and   representing test functions defined in the appropriate space,  is the domain volume, and Ā a surface where the surface traction t boundary condition is applied.We assume homogeneous Neumann boundary conditions for the momentum moments balance.Using standard finite elements for space discretization, the kinematic fields are approximated with interpolation functions.For the ease of notation, the interpolation function for each node  is formally written in a global operator   .Using this operator, the interpolated displacement component , for example, reads   =      .In this approximation,    contains the th component of the nodal displacement vector of node .The nodes are indicated with the subindex with a capital letter (⋅)  or (⋅)  .Considering that the test functions use the same interpolation operator, substituting it in the weak form Equations ( 21) and ( 22), we obtain a set of nonlinear discrete equations, with the operator  , representing the spatial derivatives of the interpolation functions and    and    denoting the residual vector components.
Because of the nonlinear material behavior, Equations ( 23) and ( 24) are treated incrementally.Equilibrium is obtained in each consecutive load increment using the Newton scheme.In each step of the Newton procedure, a linearized problem is solved for a Newton increment of the primary fields, which is collected in a vector .The matrix  is the tangent matrix, and  is a vector containing the negative residuals for each field.Equation ( 25) is solved iteratively using a preconditioned Krylov method, and the result is used to update the fields to the next Newton step.Algorithm 1 presents a pseudo-code with the overall solution process.In the algorithm, many details were omitted to keep it concise.We emphasize aspects related to the presented work, namely, the preconditioner operator setup and the usage of this operator by the linear solver.
Because of the multifield nature of the problem, the tangent matrix has a block structure of the form where the explicit expressions for the blocks in the elastic regime are given by Those expressions are required for the construction of the preconditioner in the next section.Within the above equations, the  ∈ ℝ   ×  block is related to displacement DOF and  ∈ ℝ   ×  to microrotation DOF.The values   ,   refer to the respective number of DOF.The off-diagonal blocks,  1 and  2 , are coupling blocks.The constitutive relations in this case are computed using Equations ( 14) and ( 15).
In the plastic regime, some extra terms appear due to the linearization procedure.Furthermore, the  block loses symmetry because of the nonassociated flow rule, used for modeling the behavior of cohesive-frictional materials.Also in the plastic regime, the coupling blocks do not generally satisfy the relation  1 =   2 .The explicit expressions for the blocks in the plastic regime are given by

Preconditioner development
To solve the linear system of equations (25) we adopt an iterative preconditioned GMRES solver.The preconditioner for this system is constructed based on the block factorization of the system.Following Ref.,11 a standard lower triangular block preconditioner then can be constructed as where  =  −  2  −1  1 is the Schur complement.The lower triangular preconditioner is favored instead of a block diagonal preconditioner as  is generally nonsymmetric. 10In practice, the Schur complement needs an approximation such that its inverse applied to a vector is inexpensive.This approximation is discussed in the next section.

Schur complement approximation
In the preconditioner (35), the Schur complement  requires the inverse of the  block.However, the exact Schur complement is generally dense, and its formation is prohibitive in a practical setting.Furthermore, a matrix-vector multiplication with dense factors is not desirable.Therefore, we seek an appropriate approximation for the Schur complement operator.
Because there are many possibilities to approximate the Schur complement, this step is considered application-dependent.The most common approach to deal with the Schur complement approximation is to diagonally lump the terms of the matrix to be inverted, that is,  −1 ≈ (diag ) −1 . 33,34This strategy keeps the Schur complement sparse and makes its inverse computation straightforward.On a similar note, instead of using the diagonal entries, one can compute the row or column norm and use that instead. 35 well-known method for Navier-Stokes problems is the spectral equivalence of the Schur complement with a mass matrix, which allows for an optimal approximation, 10,36 in which the eigenvalues of the preconditioned system are bound by an interval independent of the problem size.A use of a weighted mass matrix in the approximation of the Schur complement is also reported in Ref. 37 for coupled magma/mantle dynamics, and in Ref. 38 for coupled multifield problems in geomechanics.
More recently, a scheme that approximates the inverse only in a sparse pattern was presented in a general setting in Ref. 39 For our particular setting, the Schur complement approximation's accuracy increased with denser sparsity patterns, but this did not enhance the linear solver's convergence behavior despite the higher computational cost.
Another strategy, based on operator analysis for a particular problem involving thermal-fluid coupling, recognizes that  −1 takes the form of a compact perturbation of the identity matrix. 40This indicates that the approximation  ≈  can provide reasonable convergence properties and avoids the need to form the inverse of .This strategy is equivalent to a block Gauss-Seidel 41 preconditioner, which uses the lower triangular split of  to form the iteration matrix.
In the context of this work, we adopt the Schur complement approximation  ≈ .This is not, however, a generalpurpose approximation.Nonetheless, in our experiments, we found that it works well for the presented model.The approximation quality can be analyzed based on an estimate of the orders of magnitude of the blocks' entries in a linear elastic setting, Equations ( 27)- (30).From the following product between the Schur complement and the microrotation-related block inverse, one can observe that, if dense and small, the second term on the right-hand side represents a perturbation to the identity matrix.The magnitude of the entries of the displacement stiffness block  is governed by the shear modulus , see Equation ( 14), or, equivalently, by the elastic modulus .Due to the volume integral and the two shape function derivatives present in Equation ( 27), the respective entries can, furthermore, be assumed to be proportional to the finite element dimension ℎ.
The magnitude of entries in the microrotation block  is defined by two terms.The first comes from the moment of moments balance, which is the first term of the right-hand side of Equation (30).Its value is governed by the micropolar constitutive relation in Equation (15), which has the magnitude given by 2 + .To obtain these micropolar parameters in terms of , we use a relation following from the physical constraint on the material parameters, 5  .Hence,  is roughly of the same order of magnitude as .
Moreover, it is scaled in the same way as the displacement stiffness block by the element dimension ℎ.The second term on the right-hand side of Equation (30) comes from the dependence of the strain   on the microrotation   , Equation (17).In it, the alternating tensor is applied twice to the stress-strain constitutive operator, which results in a value order of   .Since no derivatives are involved, this term is scaled by the element dimension ℎ 3 , which tends to be small as the mesh is refined.
Finally, the coupling blocks have entries whose magnitude is governed by   and are scaled by ℎ 2 .Combining those results, we obtain the following rough approximation of the magnitudes of entries The notation ∼ indicates that a representative scalar entry in the matrix has the specified magnitude.Using those estimates in Equation ( 36), the entries of the second term have a magnitude of order The upper bound case   =  results in a magnitude given by 1∕(  2 d ℎ 2 + 1), which is strictly smaller than unity.When   is smaller than , this result would be even smaller.This implies that the approximation for the Schur complement  ≈  would produce an effective preconditioner.This statement is further supported by an analysis of a small-scale example problem which is documented in Appendix A. It Exemplifies a distinct clustering effect on eigenvalues of the preconditioned matrix for a critical simulation step.We also tested the common approximation using the (diag  −1 ) to form the inverse of  in the Schur complement expression, This approximation exhibits similar convergence properties as the adopted Schur complement approximation up to the peak load.However, in the softening regime, the preconditioner with the approximation in Equation ( 39) produces worse convergence behavior and requires more iterations to converge than the preconditioner with the Schur complement approximation  ≈ .The approximation in Equation ( 39) is, thus, not feasible in our context.

Preconditioner application
The action of the preconditioner in Equation (35) will require the action of the individual blocks  −1 and  −1 on general vectors.The action of the inverse to a known vector is equivalent to solving a linear system.When the linear system is approximately solved, this action is approximated.In practice, a fixed number of iterations from another iterative solver is used for this purpose.The quality of this approximation determines the effectiveness of the preconditioner.In this work, the inverse action of the blocks is approximated with a single iteration of the algebraic multigrid method (AMG) V-cycle.Hence, the preconditioner is denominated as a block preconditioner applying the AMG method for the inverse approximations of the diagonal blocks (B-AMG).Since the blocks are different, different parameters can be used during the setup of the multigrid hierarchy.
Algorithm 2 shows the procedure to apply the preconditioner operator to a vector .In the algorithm, we use the DOF from each field to split the input vector .The vector split is represented by the notation  [  ] , which indicates the vector associated with the displacement sub-block.

Numerical simulation of a compression test on sandstone based on the micropolar continuum without damage
The first application consists of the numerical simulation of a plane strain compression test on a sample of Gosford sandstone as depicted in Figure 1.In the present subsection, the displayed results are obtained using a large-deformation version of the simplified elastic-plastic model that was reviewed in the geometrically linear context in Section 2.1.Therefore, damage-induced softening is not yet included, and thus the solution can be expected to deviate significantly from the experimental results obtained for the compression test on Gosford sandstone by Ord et al. 42 Nevertheless, the results presented in this subsection deliver valuable information, as they allow us to evaluate the benefits of using the proposed block preconditioner for micropolar continua.The material parameters used are summarized in Table 1.The parameters correspond to the subset of elastic-plastic parameters that have been calibrated for the experiments in Ref. 42  micropolar damage-plasticity model, which will be used in the following Section 3. Their physical meaning was already introduced in Section 2.1.
The domain is fixed at the bottom and loaded with a uniform displacement of 4 mm at the top surface after a confinement pressure of 20 MPa is applied.Additionally, for considering the constraint imposed by the interaction between the load plate and the specimen, a rigid body constraint is approximately assumed on the top and bottom by increasing the stiffness of the material in those regions.With this constraint, the relative displacement between the load plate and the specimen is reduced.A small perturbation in one of the corners is added to trigger the direction of the shear band.
The micropolar framework was demonstrated to be effective in regularizing the instability related to the nonassociated plastic flow rule. 5Therefore, as damage-induced softening is disregarded, the employed micropolar continuum is sufficient for proper regularization in the present simplified setting.
The load-displacement curve of the structure is shown in Figure 2.There we can notice that even without damage, the model presents a slight structural softening response after the onset of localization, which can be attributed to largedisplacement effects.However, it follows from the comparison of the computed load-displacement diagrams with the experimentally determined one that consideration of the evolution of damage is necessary to reflect the actual structural behavior in the post-peak domain.
The equivalent plastic strain and microrotation fields are visualized for point C, indicated in Figures 2 and 3. From the contour plot, we can see that the plastic behavior concentrates in a narrow zone.
Concerning the discretization, quadrilateral elements with quadratic shape functions and reduced integration are used for the displacement and microrotation variables.The computation is performed for three different finite element discretizations, including a coarse, medium, and fine mesh with 7563, 29,523 and 116,643 DOF, respectively.As can be seen in Figure 2, the results are identical for all three meshes.

Method B-AMG M-AMG
Near null-space candidates   : RBM a   : Abbreviation: RBM, rigid body motions.AThree translations and three rotations for three-dimensional problems.BThree candidates as constant vector for each microrotation component.
In what follows, we demonstrate the effectiveness of the sparse Schur complement approximation during the preconditioning step of the solution process.We also illustrate the shortcomings of black-box preconditioners for multifield problems, in particular, for generalized continuum mechanical models with additional microrotation fields.For that purpose, we utilize a monolithic algebraic multigrid method as a black-box preconditioner, which is applied to the whole system.It is referred to as M-AMG in the following.Within a general AMG setup, the coarsening step uses an aggregationbased algorithm with transfer operator smoothing.The monitored parameter for the preconditioner will be the iteration count for the GMRES Krylov solver to achieve a specified relative tolerance of 10 −5 , which is used for all experiments performed.To isolate the effect of the preconditioner, the Krylov basis vectors were all kept.This means that the GMRES convergence will not degrade because of the restart option. 41The total number of iterations allowed was kept high so we can illustrate the poor performance of an inadequate preconditioner.Ideally, the iteration count should not be dependent on the discretization or a material parameter.
Algebraic multigrid methods are known as specialized tools since they have multiple user-defined options.For all numerical simulations, the same level smoother was employed for both strategies, namely one sweep of the symmetric Gauss-Seidel solver.Table 2 presents a summary of the near null space candidates used.The option to improve the near null space candidates represents a few iterations on a homogeneous system to ensure the kernel representation.We used four iterations of the Gauss-Seidel solver for all numerical simulations.The prolongator smoothing option is employed to provide convergence benefits associated with energy-minimization arguments. 44For this step, one Jacobi iteration is used.The implementation of the GMRES solver is provided by Virtanen et al. 45 and the AMG by Bell et al. 46 The material model was provided by the material library Marmot. 47he block preconditioner approach allows different AMG setups for each field.This permits using information from the field operators when choosing the near null space candidates.For the block related to displacements, the rigid body motions were used as near null space candidates. 27For the block of microrotations, the appropriate choice of near null space candidates is not that straightforward.Since the sub-block D is nonsingular, no obvious candidates exist in this case.The vectors contained in this near null space thus have to be interpreted as slow converging error modes rather than as rigid body motions.Our choice is a constant vector of one's for each micro rotation component.It was motivated by an investigation of the eigenvectors of the sub-block of microrotations of an element stiffness matrix.We found that the mode associated with the smallest eigenvalue was approximately a constant vector.For the monolithic preconditioner, a constant vector is used, since there is no obvious near null space candidate when multiple fields are considered simultaneously.The number of iterations for each GMRES call is accumulated during a time step.This value is normalized by the number of nonlinear iterations in the time step and plotted for each step in Figure 4. Figure 4B shows the iteration count when the tailored block preconditioner developed is used for the computation with the assumed ratio   ∕ = 0.1.The same curves are shown in Figure 4A together with the results when the black-box monolithic algebraic multigrid method preconditioner is used.From this figure, we can see a high sensitivity to the problem size, which evidences the problem with this type of preconditioner in the context of multifield models.Not only the number of iterations is higher, but the total computing times are of order 10 times larger 3 .Another remark can be made on the computed structural response: as it approaches the peak load, the system becomes increasingly ill-conditioned, represented by the increase in the number of iterations required.

EXTENSIONS FOR MICROPOLAR CONTINUUM MODELS ENHANCED BY A NONLOCAL DAMAGE FORMULATION
In this section, the developments from the previous section are extended to models accounting for material degradation by using an additional nonlocal damage field.Hence, the resulting systems are commonly larger in size, and feature an extended block structure.Furthermore, the developed numerical methods now have to deal with more demanding physical processes, as the extended model exhibits a significantly more pronounced softening response on the structural level, including the formation of thinner, more localized failure zones.The change in the system properties needs to be considered for the development and analysis of the block preconditioner strategy.
We begin with a brief presentation of the model formulation related to damage.Then, the discrete treatment with finite elements is discussed.Subsequently, an extension of the 2 × 2 block preconditioning strategy, developed in the last section, to 3 × 3 block systems will be described.This approach for including the nonlocal damage block has certain similarities to approaches established for phase field models, however, the context here is wider due to the presence of the micropolar model.The proposed approach is analyzed for the model problem already used in the previous section.However, due to the extension of the model by including the evolution of damage, the model is able to better represent the response of the actual experiment.

Nonlocal damage formulation
In the presently considered case of a micropolar isotropic damage-plasticity model, degradation processes of the material are accounted for by a single scalar damage parameter , which is used for representing both, the degradation of the stress as well as the degradation of the couple-stress.This damage parameter can be used to convert nominal (couple) stresses to effective (couple) stresses.Nominal stresses are to be interpreted as forces per area, while effective stresses can be understood as forces per intact area of the material.The damage parameter is assumed to depend on the nonlocal damage driving field αd , according to the constitutive relation This corresponds to the assumption of an exponential softening law, with  f denoting the softening modulus.
The nonlocal character of damage is modeled with the screened Poisson equation presented in Equation (11).When modeling failure in quasi-brittle materials, nonlocal damage formulations are known to provide regularized results, that is, results independent of the employed discretization.The local damage driving variable in Equation ( 11) is driven by the accumulation of the inelastic volumetric deformation.The latter can be obtained as part of the plasticity model, which is formulated in terms of effective stresses.In the special case of a geometrically linear context, the evolution of this variable is governed by the rate of the volumetric plastic strain weighted by a function accounting for different behavior in tension and compression.In this work, we are interested in the final algebraic equation from the discretization of Equation ( 11), thus the details of the constitutive model are omitted and we refer to Ref.5 for detailed information.

The discrete system
By analogy to Section 2.2, the discrete equation related to damage is obtained with the finite element method.Once more, for the sake of brevity, the following presentation is restricted to the geometrically linear context.Following the same steps as in the previous section, we start by writing the partial differential equation (11) with  αd as a test function defined in an appropriate space.Using the same interpolation operator as for the displacement and microrotation, we obtain the discrete version, Furthermore, the residual vector components for the balances of linear momentum and momentum moments, ( 23) and (24), can be formed also for the micropolar nonlocal damage continuum.However, stresses and couple stresses in these equations are now to be interpreted as nominal stresses and couple stresses, respectively.Thus, they depend also on the damage parameter.When these three residual equations are considered in the linearization procedure, the tangent matrix will have additional blocks with explicit expressions for the latter given by The remaining four blocks ,  1 ,  2 , and  are formally given by the same equations as for the micropolar model, that is, by Equations ( 31)- (34) in the geometrically linear context.However, as in the definition of the residuals, the stress and couple stress derivatives in these equations are now derivatives of nominal stress and nominal couple stress, respectively.A distinct feature of this system is that it undergoes significant changes in the material behavior (in addition to plastic deformation) during the simulation process.Coupling blocks related to damage become nonzero when the damage process initiates.

Preconditioner development
We will extend the preconditioner, developed in the previous section for the micropolar continuum only, to introduce the new nonlocal damage field.A lower-triangular block-preconditioner for the matrix in Equation ( 43) can be constructed the same way as it was done for the two-field system.Assuming the matrix  is invertible, block-factorization of it yields a full-triangular preconditioner  −1  , The matrix  1 is the first-level Schur complement, and  2 represents a short notation for The term  2 is a second-level Schur complement, which is given by The full-triangular preconditioner in Equation ( 49), when applied exactly, would result in a preconditioned matrix whose spectrum consists of the unit eigenvalue only.However, the preconditioner  −1  is not suitable because it demands operators not easily available, making it not practical.We need, therefore, algebraic approximation for  −1 ,  −1 1 , and  −1 2 to reduce the computational cost of preconditioning the system.There are many ways to approximate those operators, we investigated three variants.The first strategy under consideration is to use the first-level Schur complement  −1 1 ≈  −1 to construct a sparse second-level Schur complement by replacing the inverse operators with a diagonally lumped one.In this situation, the second-level Schur complement, from Equation (52), would assume the form The second strategy is to further simplify this approximation of the second-level Schur complement such that  2 ≈ .Assuming that  −1 2 ≈  −1 , and in conjunction with the suggested approximation for the first-level Schur complement  −1 1 ≈  −1 , we can formulate a preconditioner M−1  as follows: In this form, a preconditioner application to a vector would need four block inverses approximations.Considering the application of this preconditioner to a known vector , we want the result The first approximation of the diagonal blocks inverses comes from the first block-row of equations.This means we are solving In which the known vector  1 is a split of  corresponding to the first field and  1 a split of the solution vector we want.This system is approximately solved using an iterative solver, for example, algebraic multigrid method.The second inverse approximation arises in the second block-row of equations, where we need to compute with  2 =  2 −  2  1 as an intermediary value.The third block-row of equations yields a third inverse approximation, where  3 =  3 −  2  1 −  2  2 +  2  4 as another intermediary value.In  3 , one needs the fourth inverse approximation to compute Alternatively, if we inspect the system before damage initiates, we notice that the coupling blocks  2 and  2 have only zero entries.Using this fact, we can construct the left lower-triangular preconditioner  −1 based on the system before damage initiates, After damage initiates, the application of this preconditioner would not generate a spectrum consisting of the unit eigenvalue only, as illustrated in Appendix A, because of the additional coupling with the nonzero blocks  2 and  2 .However, the preconditioner is cheaper to apply than the previously defined versions.
Although the use of a full-triangular preconditioner is often recommended for nonsymmetric problems, 10 we found based on numerical experiments that for the problem considered in the present publication, the associated gain in iteration count using the full triangular preconditioner in Equation ( 54) is not worth the extra computational effort in comparison with the variant in Equation (59).The same can be said when we use the second-level Schur complement approximation in Equation (53) with the full triangular preconditioner from Equation (49) and the first-level Schur complement approximation  −1 1 ≈  −1 .It also produces similar convergence properties with additional computational costs.Therefore, we use the simpler preconditioner shown in Equation (59) in all of the following computations.

Numerical simulation of a compression test on sandstone based on the micropolar continuum model enhanced by a nonlocal damage formulation
This case study deals with the numerical simulation of the compression test on Gosford sandstone performed in Ref. 42  It is an extension of the previous numerical study in Section

F I G U R E 5
Reaction force as a function of the vertical applied displacement at the top surface of the domain for different mesh sizes for the ratio   ∕ = 0.1 (left) and for the fine mesh assuming different ratios   ∕ (right).
damage is present.The extended model comes along with two additional material parameters, the softening modulus   and the length scale parameter for damage  d , for which values are provided in Table 3.The values of the other parameters, shown in Table 3, are equivalent to the parameters used in the numerical study in Section 2.4.For the numerical simulation, a two-dimensional model is presented first.Further results for the three-dimensional setting, for which the microrotation field has three components, will be included in the next section.The geometry, boundary conditions, and material parameters of the compression test are the same as in Section 2. 4 and specific details about it can be found in Ref. 32 In regards to the discretization, the quadrilateral elements from the previous simulation are extended by including the additional nonlocal damage-driving field.The linear solver preconditioning strategy is tested for three different discretizations, denoted as coarse, medium, and fine with 10,084, 39,364, and 155,524 DOF, respectively.
The response of the structure is shown in Figure 5, where the reaction force is plotted as a function of the applied vertical displacement at the top surface.Note that as the model discretization is refined, the solution converges to a unique response, which corresponds very well with the experimentally determined one.This stable result is attributed to the localization limiters employed.The effect of the material parameter ratio   ∕ is examined as well in Figure 5. Therein we observe that for the fine mesh, the load-displacement curves are similar for the different ratios of the coupling modulus over the shear modulus.We also notice that for a small   ∕ ratio, the response is slightly more brittle, which may be associated with the reduced effect of the microrotations as a localization limiter.The post-peak softening response is represented by the nonlocal damage formulation.
Also in Figure 5, two instants were selected as represented by arrows labeled A and B. The damage field for those instants is plotted on the deformed structural shape in Figure 6.There we notice the localization of the damage process right after the peak load.Figure 7 shows the plastic strain, microrotation and nonlocal damage fields at the last step, which is indicated with the label C in Figure 5.There we can also see the deformed shape of the structure and the shear band.

Evaluation of preconditioner performance
The applied displacement is subdivided into 500 steps.Each step can involve multiple nonlinear iterations.Each nonlinear iteration is solved with the GMRES method, as in Section 2.4. Figure 8B shows the accumulated number of linear solver iterations normalized by the number of nonlinear iterations for each displacement step using the developed three-field B-AMG.One can see the effect of the damage process on the effectiveness of the preconditioner.As the material further degrades, the system experiences changes which impair the performance of the preconditioner, facilitated by the fact that it was constructed based on the system before damage initiates.Figure 8A presents a comparison between the monolithic (A) (B) Comparison of GMRES iterations for different material parameter ratios   ∕ using the proposed block preconditioner B-AMG for the medium mesh (left) and for the fine mesh (right).

TA B L E 4
Average number of linear solver iterations for the whole simulation steps.strategy M-AMG and the block strategy B-AMG.We can see a considerably higher iteration count and larger oscillation in the monolithic strategy, which evidences the advantages of the proposed block preconditioner for the problem.

B-AMG
The simulation is repeated for different discretizations and different material parameter ratios   ∕. Figure 9 shows a subset of the results.Again, the accumulated number of linear solver iterations normalized by nonlinear iterations is plotted.In Figure 9a, the results obtained for three different ratios of the coupling modulus to shear modulus   ∕ are shown for the medium discretization.In Figure 9B, the respective results are shown for the case of the fine mesh.The material parameter ratio   ∕ has a small effect on the linear solver performance, except at the onset of localization.At the beginning of localization, with a smaller coupling modulus, the GMRES method requires more iterations to converge.This effect increases as the mesh is refined.
Table 4 shows the average number of linear solver iterations for different mesh sizes and different material parameter ratios   ∕ for the proposed preconditioner and for the black-box alternative.Note that as the number of DOF increases the average number of iterations required also increases for both preconditioners.However, the sensitivity is considerably smaller for the proposed preconditioner.

A PROBLEM-ADAPTED SETUP STRATEGY FOR THE DEVELOPED THREE-FIELD BLOCK PRECONDITIONER
The computed structural response, depicted in Figure 5, has a very characteristic pattern, which features a distinct peak load followed by rapid structural softening.This distinctive response can be exploited in the solution procedure.The present section presents an optimization approach to exploit this feature.It will be evaluated in 2D and 3D numerical numerical simulations for the compression test considered before.

Development of a problem-adapted update strategy of the preconditioner
In this subsection, a new approach for a selective update of the preconditioner is presented.The idea aims at optimizing the preconditioning setup step, more specifically, the setup of the algebraic multigrid hierarchy.The setup time cost can be a significant part of the total linear solver time, especially in simulations that consist of multiple sequential steps.
F I G U R E 1 0 Scheme of the problem-adapted strategy for the three-field preconditioner.
The problem-adapted strategy is developed on two key observations.The first one is related to those parts of simulations in which the system only slightly changes between consecutive steps.This allows exploring the reuse of the algebraic multigrid hierarchy computed for the system in previous steps.The second observation is related to the fact that even if the system changes between consecutive steps, the changes during the nonlinear iterations are small.
Considering those observations, one can conclude that in the prepeak domain of the load-displacement curve, the AMG setup can be reused for subsequent preconditioners without impairing the linear solver performance.At peak load, the onset of localization of strains and damage occurs and, as a consequence, the system becomes ill-conditioned, which is reflected in the iteration counts of the linear solver.Thus, in the post-peak domain, reusing the AMG setup from the beginning of the simulation is not effective anymore as it slowly deteriorates the linear solver performance.This performance decline is related to the damage evolution.
Since we are dealing with a nonlinear problem, multiple nonlinear iterations are required for each time step.From the second observation, reusing the AMG hierarchy within a step does not impair the linear solver performance in terms of the number of iterations required.Thus, there is no need to recompute the hierarchy at each nonlinear iteration.
Figure 10 shows a schematic representation of the strategy developed.

Evaluation of the problem-adapted setup strategy
The new problem-adapted strategy  B-AMG for the developed block preconditioner (with  indicating the problemadapted strategy) is evaluated for the simulation of the compression test in Section 3.4 by comparing the time spent during the setup of the preconditioner for the two methods.The usual adoption of a custom block preconditioner requires its setup based on the system matrix.

Evaluation of the problem-adapted update strategy for 2D simulations
For the present case, the preconditioner setup is performed every time the linear solver is called, which happens at every step of the simulation and every nonlinear iteration.The problem-adapted strategy exploits this frequent setup operation to reduce computation time.The strategy divides the simulation into before and after the peak load.Before the peak load, the preconditioner setup is reused on the subsequent time steps.After the peak load, the AMG setup is computed only at the first nonlinear iteration.
In Table 5, the proportion of setup and linear solver wall times are presented.The same result can be viewed graphically in Figure 11.From the reported measures, the time spent in the setup for the usual approach is in the range of 10%-25% of the total time spent solving the linear system.We observe a reduction in the setup time of order of 3 between the usual application of a block preconditioner and the problem-adapted strategy.This reduction in setup time has no negative impact on the convergence properties.Therefore, the time savings from not recomputing the AMG hierarchy reduce the total computation time.However, the benefit of the problem-adapted strategy is somewhat diminished when a smaller  ratio   ∕ is used.This is due to difficulties in solving the system at the onset of localization, indicated in Figure 9 by a pronounced peak in the iterations count.

Evaluation of the problem-adapted update strategy for 3D-simulations
The three-dimensional numerical simulation of the compression test on Gosford sandstone represents an extension of the plane strain simulation.The main difference from the plane strain case is that the microrotation field now has three components.This has an impact on the block preconditioner strategy for the approximation of the inverse of the microrotation block.The domain and boundary conditions of the 3D model are represented in Figure 12.The model is free to move on the  −  planes.Additionally, one corner node is completely fixed to prevent a rigid body translation.A confining pressure of 20 MPa is applied on the surfaces in the  − -planes.A vertical displacement of 4 mm is specified on the top surface.A small perturbation is applied on a top corner to trigger the shear band direction.For layers of 10-mm thickness on the  top and bottom of the specimen, the material stiffness is increased to represent the constraint between the load platens and the specimen.The domain is discretized with 20 nodes quadrilateral elements with reduced integration.Two mesh resolutions are used: a coarse and a fine mesh with 70,301 and 491,127 DOF, respectively.The computed structural response is similar to the two-dimensional case shown in Figure 5.
The objective is to check if the performance of the problem-adapted preconditioning strategy in a larger threedimensional problem is similar to the smaller two-dimensional case.The numerical simulation is performed with the problem-adapted scheme developed,  B-AMG, and the regular procedure with recomputed preconditioner setup at every nonlinear iteration, B-AMG, is used as a reference.
The experiment is performed for three different ratios   ∕ and the wall time of the execution is monitored.The proportion of the preconditioner setup wall time to linear solver wall time and the total linear solver time are summarized in Table 6 and are displayed in Figure 13.The results demonstrate that the setup time can be reduced by a factor of 3 when comparing a regular preconditioner approach and the developed problem-adapted strategy.This reduction occurs without impacting the solver convergence properties which results in overall time savings.However, this gain is not noticeable when a smaller coupling modulus is used, as shown with the ratio   ∕ = 0.01.Which is the same effect experienced in the two-dimensional case.
In Figure 14, we present the accumulated iterations count from the linear solver, normalized by the nonlinear iterations, for the problem involving a material parameter of   = 0.1.Following this, we examine the sensitivity of the linear solver to the material parameter.We perform this analysis separately for the coarse mesh (Figure 15A) and the fine mesh Comparison of GMRES iterations for different material parameter ratios   ∕ using the proposed block preconditioner B-AMG for the coarse mesh (left) and for the fine mesh (right).
(Figure 15B).Similar behavior was observed in the 2D case as shown in Figure 9, exhibiting some sensitivity to the mesh size as the simulation progresses within the softening regime.

CONCLUSION AND OUTLOOK
Regularization of numerical models aiming at the simulation of localized failure phenomena is crucial to ensure meshindependent results.Such regularization can be achieved using generalized continuum approaches, which introduces additional unknowns and additional balance equations.For localized failure, the resulting system of equations is highly nonlinear and in most cases, it is characterized by a large number of DOF.Its is commonly computed by solving a series of large and sparse linear systems.It is common practice to perform a Newton-Krylov approach, that is, to obtain these linear solutions using a Krylov-subspace-based iterative method.Therefore, the simulation of localized failure phenomena heavily relies on the efficiency and robustness of the Krylov iterative solvers, for which a preconditioning step ensures and accelerates convergence, in the optimal case, independently of the discretization.For the generalized continuum models considered in the present publication, the linear systems exhibit a physics-based block structure, which was exploited for the preconditioner construction.To this end, preconditioning strategies were developed for multifield problems characterized by a micropolar continuum formulation and they were further enhanced by a gradient-enhanced nonlocal damage model.
In a first step, a preconditioner for linear systems resulting from problems simulated with micropolar models, which involve only displacements and microrotations, was constructed based on a block system factorization of the corresponding 2 × 2 block system.It employed a sparse approximation of the Schur complement based on the microrotation block, which was deduced from an analysis of the sub-blocks' orders of magnitude.For the application of the resulting triangular block preconditioner, the inverse approximations for the diagonal blocks were computed by means of the algebraic multigrid method.
The appropriate choice of near-null space candidates is crucial for the effectiveness of algebraic multigrid preconditioners.Therefore, a first benefit of the block preconditioning approach was found to be that it was more straightforward to select the appropriate null space candidates for the individual diagonal blocks, then it was for the coupled system, for which null space candidates can include combinations of displacement-based macrorotations and microrotations.More specifically, for the displacement sub-block in the block preconditioning approach, the rigid body modes for displacementbased elasticity problems are used as near null space candidates.For the microrotation block, the choice of near null space candidates as constant microrotations worked well.
Results of numerical tests on a micropolar problem indicate that the proposed preconditioner ensures rapid convergence in the considered test case.Iteration counts show only a mild dependency on the mesh resolution, and the block-preconditioned approach significantly outperforms the black-box monolithic method.
In a second step, the previously developed preconditioner was extended for micropolar continua enhanced by a nonlocal damage formulation.This extension was based on the structure of the 3 × 3 block system before damage initiation, and can be understood as a block diagonal extension of the previously developed 2 × 2 preconditioner with respect to the additional nonlocal damage field.
Again, all diagonal block inverses were approximated by means of the algebraic multigrid method, this time using the usual constant vector, in which all entries are equal to one, as near null space candidate for the additional nonlocal damage field.
The preconditioning strategy for the three-field formulation was tested for a Gosford sandstone compression test, an experiment featuring localized failure in a narrow shear band accompanied by large deformations.The computations again showed that the block-preconditioned approach significantly outperforms the black-box monolithic method.However, for the three-field case, an increase in the iteration count can be observed for the proposed block preconditioning approach as failure progresses.This indicates the limitations of the present approach, which all couplings between nonlocal damage and the other two fields are neglected in the preconditioner.This choice was motivated by the requirement of a practical preconditioner.By using the system before damage as the basis to construct the preconditioner, we avoid an expensive second-level Schur complement construction.
In a third step, a novel update strategy was implemented by combining aspects from the simulation stage and nonlinear iteration into the linear solver.The algebraic multigrid setup is selectively updated during the solution process.Before the simulation reaches the peak load, the AMG setup is reused on the subsequent preconditioner calls.This reduces total computation time by avoiding unnecessary preconditioner setup.The nonlinear aspect is incorporated by only computing the AMG setup on the first nonlinear iteration.
Future research topics may include the parallel computation of the developed strategy.Another potential area for investigation is the interaction between the tolerances of the nonlinear and linear solvers.Instead of using the same small tolerance throughout the whole simulation, an adaptive tolerance may provide performance benefits in the context of inexact Newton methods. 33

A C K N O W L E D G M E N T S
This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 847476, which is gratefully acknowledged.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.In this appendix, we conduct an eigenvalue analysis of the preconditioned matrix for a small example problem.The eigenvalues are computed for the system matrix resulting from a specific increment of the two-dimensional coarse problem in Section 3.4 involving three fields: displacements, microrotations and nonlocal damage.For illustration purposes, we selected the most critical increment, which is at the onset of localization.Preconditioning is done using the block preconditioner presented in Equation (59), which is repeated here for convenience,

O R C I D
The preconditioned system  −1  is obtained by explicitly forming the inverse of the sub-blocks ,  1 , and .First, the analysis is performed using the exact expression for the Schur complement, that is,  1 =  −  2  −1  1 .The resulting eigenvalues of the preconditioned system are plotted for three different simulations considering distinct ratios   ∕ in Figure A.1.Also in the figure's legend, we indicate the ratio between the maximum and minimum positive real eigenvalues and the absolute value of the smallest real eigenvalue.Since we are using the exact expression for the Schur complement, and for the inverses, this clustering represents an "ideal" scenario.When applied to the system before damage initiates, the spectrum consists of the unity eigenvalue only.In practice, this "ideal" clustering can not be achieved since the Schur complement and the inverses of the sub-blocks are approximated.From the distribution of the eigenvalues, we can notice that as the ratio between the coupling modulus and shear modulus increases, the smallest real eigenvalue gets closer to zero.Second, if we now use the block-preconditioner from Equation (59) with the Schur complement approximation  1 = , we get the eigenvalues distribution shown in Figure A.2. From the figure, we can observe that this Schur complement approximation also produces a clustered spectrum.Here the behavior related to the material property is also similar, as the coupling gets larger, the smallest real eigenvalue gets closer to zero.

F I G R E A . 1
Eigenvalues of the preconditioned system using the exact Schur complement expression for different coupling modulus ratios   ∕.

1
Simulation overall procedure with Newton-Krylov solver procedure run simulation while  ≤   do Increment the load while not in equilibrium do Assemble linear system  =  Setup preconditioner operator  −1 Call linear solver GMRES(,  −1 , ) ▹ Application of  −1 is shown in Algorithm 2 Update the solution and check equilibrium end while Increment the time  =  + Δ end while end procedure t = 2 b and  t = √   , with  t as a characteristic length for torsion.Combining those, we find  = 4 2 b .The parameter  is obtained also with a relation following from the physical constraint using the full F I G U R E 1 Domain and boundary conditions for the two-dimensional problem.F I G U R E 2 Reaction force as a function of the vertical displacement at the top surface of the domain for different mesh sizes and the material parameter ratio   ∕ = 0.1.

F I G U R E 3
Plastic strain (left) and microrotation field (right) displayed on the deformed structure at the last time step, which is indicated by the label C with an arrow in the load-displacement curve in Figure2.TA B L E 2 Summary of near null space candidates for the AMG.

4
Comparison of the performance of the monolithic black-box preconditioner M-AMG with the proposed block preconditioner B-AMG for the micropolar continuum model with the ratio   ∕ = 0.1 in regards to the accumulated linear solver iterations related to the number of nonlinear iterations per time step (p.t.s) for the coarse, medium, and fine discretization (A); the results for the proposed block preconditioner B-AMG are plotted again on a different scale (B).

F I G U R E 6 5 .F I G U R E 8
Distribution of the damage parameter  for two different instants, A and B, which are indicated in the load-displacement curve in Figure 5. F I G U R E 7 Plastic strain (left), microrotation (middle), and nonlocal damage (right) fields for the displacement, indicated by the label C in Figure Comparison of the performance of the monolithic black-box preconditioner M-AMG with the proposed block preconditioner B-AMG for the model based on the micropolar continuum (  ∕ = 0.1) with nonlocal damage in regards linear solver iterations count per time step (p.t.s) for the coarse, medium, and fine discretization (A); the results for the proposed preconditioner B-AMG are plotted again on a different scale in (B).

F I G U R E 1 2
Domain and boundary conditions for the 3D numerical model.

F I G U R E 1 4
Sensitivity of the preconditioner to the discretization for the model with   ∕ = 0.1.
Application of the preconditioner to a vector procedure apply    ,   ← getDofIds() ▹ Field-split of the system  [  ] ← AMG(,  [  ] )  [  ] ← AMG(,  [  ] −  2  [  ]) Material parameters used in the numerical simulation of a plane strain compression test, based on the elasto-plastic micropolar formulation without damage.different∕ ratios are assumed: 0.01, 0.1, and 1, which are in the range of values reported in Ref.43 A L G O R I T H M 2 aThree in the weak form, Preconditioner setup wall time proportion to linear solver wall time (%) and total linear solver time (minutes) in the parentheses.Time spent on setup for AMG hierarchy for the 2D simulation.
Preconditioner setup wall time proportion to linear solver wall time (%) and total linear solver time (hours) in the parentheses.Time spent on setup for AMG hierarchy for the 3D simulation.