On non‐stationary polarization methods in FFT‐based computational micromechanics

Polarization‐type methods are among the fastest solution methods for FFT‐based computational micromechanics. However, their performance depends critically on the choice of the reference material. Only for finitely contrasted materials, optimum‐selection strategies are known. This work focuses on adaptive strategies for choosing the reference material, details their efficient implementation, and investigates the computational performance. The case of porous materials is explicitly included. As a byproduct, we introduce a suitable convergence criterion that permits a fair comparison to strain‐based FFT solvers and Eyre–Milton type implementations.


State of the art
As a typical model for computational homogenization 1,2 on a periodic cell Q and at small strains, we are concerned with the optimization problem seeking the total strain field on the cell Q minimizing the average elastic energy W( ) = ⟨w(⋅, )⟩ Q under the kinematic constraints  which encode compatibility and fixed macroscopic strain . A prototype for primal methods for solving the problem (1) is projected gradient descent, which iterates for step sizes s k , where ∇W denotes the gradient of W and   refers to the orthogonal projection onto the kinematic constraint . More explicitly, for the L 2 inner product, we are led to the basic scheme 3,4

of Moulinec and Suquet
where is the average strain, C k is an elastic reference material and the Γ k -operator is defined by Γ k = ∇ s (div C k ∇ s ) −1 div in terms of the symmetrized gradient ∇ s and the associated divergence div. We denote by = w∕ the stress operator, which maps a strain field to the associated stress field. The basic scheme (2) is the prototype of a primal method, that is, a method where the kinematic constraints are satisfied at every iteration. Extensions of the basic scheme lead to fast and robust solution methods for FFT-based computational micromechanics. These include conjugate gradient methods, [5][6][7] fast gradient methods 8,9 as well as Newton, [10][11][12] and Quasi-Newton methods. [13][14][15] Please note that purely dual methods may be developed as well. 16,17 These are, however, conceptually similar to the primal methods.
An alternative to the primal methods are primal-dual methods, where primal and dual feasibility, that is, kinematic compatibility and stress equilibrium, are satisfied only at convergence, in general. Michel et al. 18,19 suggested to rewrite the optimization problem (1) in split form where  denotes the indicator function of the kinematic constraint, attaining the value 0 if ∈  and +∞ otherwise. For fixed reference material C 0 , the associated augmented Lagrangian function reads where is the field of Lagrange multipliers. Michel et al. 19(sec.4.1.2) suggested to utilize the associated alternating direction method of multipliers (ADMM) 20-22 , Upon convergence, the two strain fields coincide, that is, = e holds, and the field represents the associated stress field = ( ). Please note that the second line of ADMM (4) involves an implicit stress update, that is, the mapping e  → (e) + C 0 ∶ e needs to be inverted. This contrasts with the basic scheme (2) and its descendants which involve an explicit stress update  → ( ).
We distinguish the stress operator and the stress field, which we denote by . Upon convergence of the algorithm (4), the stress field arises from the strain field via the stress operator = ( ). Taking care of this distinction in the notation helps to clarify the structure of polarization-type methods, as they crucially rely upon implicit stress updates.
For linear elastic material behavior, Monchiet and Bonnet 23 introduced a family of methods iterating on the stress polarization P = ( ) + C 0 ∶ . As pointed out by Moulinec and Silva, 24 these methods included the Eyre-Milton method 19,25 and, for the linear elastic case, the ADMM (4) as special cases. 23,24 For linear elastic material behavior and a finite contrast between the phases, Moulinec and Silva 24 showed convergence of polarization schemes for any reference material (and reasonable damping coefficient). Also, in the class of reference materials proportional to the identity, C 0 = 0 Id, they identified the optimum convergence rate at 0 = √ − + and vanishing damping, where ± refer to the minimum and maximum eigenvalues of the elastic stiffnesses of the phases.
Building upon earlier work for nonlinear polarization methods, 26,27 Schneider et al. 28 considered nonlinear polarization methods. Based on optimum parameter selection for the Douglas-Rachford splitting 29 (equivalent to ADMM for our case), the parameter choice 0 = √ − + could also be established for nonlinear material behavior, where ± refers to the modulus of strong convexity and the Lipschitz constant of the stress operator . Furthermore, Schneider et al. 28 generalized previous work 19,27,30 and showed that the implicit stress update (4) can be evaluated explicitly for small-strain inelastic materials where the stress is related to the elastic strain be Hooke's law.
Renewed interest in polarization methods arose due to possible advantages in parallelization 31 and for treating nondifferentiable energies. 32,33 Still, treating microstructures with pores remained out of scope for polarization methods. For a start, it was not even clear whether polarization methods converged for porous microstructures, or, to be more precise, for which reference materials convergence could be ensured. Secondly, the simple rule 0 = √ − + for choosing the reference material does not make sense if − = 0.
For primal methods, it was realized that the original spectral discretization of Moulinec and Suquet 3,4,34 was numerically unstable for porous microstructures. By changing the discretization to finite difference 35,36 or finite element schemes, 37,38 convergent computational methods could be obtained. In particular, the convergence speed of the solution methods were reduced to judicious parameter selection. 7, 39 Anderson acceleration 40,41 is a Quasi-Newton method applicable to generic fixed point schemes that was successfully applied to primal FFT-based methods [13][14][15] . Combined with polarization schemes, the resulting Anderson-accelerated polarization schemes 42 showed robust and fast convergence, also for porous microstructures. Unfortunately, Anderson acceleration comes with a severe memory demand. For this reason, renewed interest was sparked concerning an adaptive selection of the reference material in polarization methods.

Contributions
In recent years, primal FFT-based computational homogenization methods have experienced renewed interest. Solution strategies are sought which unite fast convergence with a low memory footprint. For this reason, modifications and extensions of the basic scheme (2) of Moulinec and Suquet 3,4 were studied. [5][6][7][8]39 One of the insights of these works is that it is generally advantageous to select the solver parameters adaptively instead of requiring a clever user input. For instance, the Barzilai-Borwein method 39 is nothing but the basic scheme where the reference material is updated for every iteration based on the Barzilai-Borwein rule. 43 To the best of the author's knowledge, adaptively changing the reference material has not been considered for polarization methods, and the work at hand is aimed at closing this gap. We consider the general iterative scheme involving a damping (or relaxation) factor ∈ [0, 1) and a reference material C k that may change from one iteration to the next, as for the basic scheme (2). We consider mixed "boundary conditions" 44 and D k = (C k ) −1 . More precisely, the mixed boundary conditions are given in terms of a pair (P, Q) of complementary and orthogonal projectors and prescribed partial strains and stresses which satisfy the compatibility conditions The algorithm (5) represents a (more or less trivial) extension of the one considered by Moulinec and Silva 24(eq.13) to nonlinear material behavior and mixed boundary conditions, with the parameter conversion = 2(1 − ) and = . By a similar reasoning as in Moulinec and Silva, 24 it is not difficult to see that the algorithm (5) reduces to the classical ADMM scheme (4) for = 1 2 and C k ≡ C 0 . The damping coefficient ∈ [0, 1) allows to represent, both, the Eyre-Milton method 25 and the Monchiet-Bonnet scheme 23 as special cases, see Monchiet and Bonnet 23 and Moulinec and Silva 24 . More generally, the abstract nonlinear polarization schemes considered in Schneider et al. 28 are included. The reason for considering the form (5) rather than the polarization form 28 is that the former allows for adaptively changing reference material in a more natural way. Indeed, for fixed strain field * , but changing reference material C k , the polarization field P k = ( * ) + C k ∶ * changes, as well.
How to choose the reference material C k in the polarization scheme (5) is treated in Section 2. We discuss theoretical convergence results, available in the literature, specialized to our case. We also present general strategies for selecting the reference material in ADMM, including the classical residual-balancing strategy, 45 a procedure 46 developed for the non-stationary Douglas-Rachford scheme, and a Barzilai-Borwein stepsize. 47,48 In Section 3, we discuss the efficient implementation of the non-stationary polarization method (5), and also provide details on evaluating the residual during the iterations. In particular, comparability to strain-based solvers and the polarization schemes 28 is ensured. Section 4 contains computational examples, evaluating the benefits of non-stationary polarization methods, but also demonstrating their limitations.

Convergence statements
In this section, we provide a guide to convergence assertions for the ADMM (5) with adaptive relaxation parameter k ∈ [0, 1) and a reference material C k = k Id proportional to the identity, that is, a linear elastic reference material with vanishing Poisson's ratio and shear modulus k = k ∕2. The results may be extended to truly isotropic reference materials provided the phases are isotropic as well. For the work at hand, we restrict to the choice C k = k Id. We assume that the stress operator arises from a (possibly condensed) free energy density w which we assume convex and continuously differentiable in the strain variable, that is, is assumed to hold for any square-integrable strain field. Furthermore, we suppose that w satisfies a growth condition, 49(sec.2.1) s.t. the optimization problem (1) admits a solution.
For the general case, suppose that the relaxation parameters k and the coefficients k satisfy the conditions: 1. The relaxation parameter is restricted as follows: 0 < k ≤ 1∕2 for all k.

The summability condition
Then, for any initial point, the ADMM algorithm (5) converges in the sense that k − e k goes to zero in L 2 and the effective free energy at e k converges to a minimum. We refer to Xu et al. 48(sec.5) for a proof. Let us briefly comment on the conditions. The first condition excludes the case k = 0, that is, damping is necessary for the general case. Condition 2 ensures that the shear moduli of the reference materials are bounded. Interpreted in terms of a step size, 28(sec.3) this ensures that the step size does not go to infinity. Condition 3 ensures that decreases in the reference constant k have a sufficiently fast decay. There is a version for increasing k , 48(Ass.2) but this is not relevant for our problem. Indeed, the linear elastic stiffness typically provides a (pessimistic) upper bound for the reference material.
Please note that, for this general case, only a slow convergence rate can be expected. Indeed, the considered scenario covers also degenerate cases such as porous materials with fractal interfaces. For the general case, and fixed parameters, ADMM converges with a 1∕k-rate. 50,51 This rate is similar to what may be achieved by projected gradient descent. 49 After treating the general case, let us consider the best-case scenario. 28,29 More precisely, suppose that there are positive constants − and + , s.t. the stress operator is uniformly − -strongly convex and + -Lipschitz-continuous for the Frobenius inner product || || = √ ∶ on Sym(d). Then, for any reference material C k ≡ C 0 = 0 Id with 0 > 0 and relaxation parameter k ≡ ∈ [0, 1), the iterates ( k , e k , k ) produced by the ADMM algorithm (5) converge to the (unique) fixed point ( * , * , ( * )). Moreover, the estimate holds for the stress polarizations P k = k + 0 e k and P * = ( * ) + 0 * .
The fastest convergence rate is achieved for = 0 and 0 = √ − + . In particular, for fixed precision, the number of necessary iterations is proportional to the square root of the material contrast √ + ∕ − . The linear elastic case was discussed by Moulinec and Silva. 24 Please keep in mind that the choice = 0 and 0 = √ − + leads to the fastest overall convergence rate. This does not necessarily correspond to the minimum iterations for prescribed tolerance. Rather, the rationale behind optimizing the convergence rate is that, for high accuracy, the number of required iterations depends on the convergence rate, in the first place.
The contraction estimate (8) was established for the Douglas-Rachford rewriting 52 of the ADMM, which is more convenient to analyze. The assertion (8) reveals a few things. For a start, zero damping = 0, is permitted, and even favorable. This possibility was excluded for the general case treated earlier. In computational practice, however, small positive values of are typically preferred over = 0.
Secondly, by the same reasoning as for obtaining the estimate (8), it is possible to derive the bound Thus, as long as the sequence k is convergent to * ∈ [0, 1), the adaptive ADMM converges linearly, as well. The third insight from the estimate (8) is that a clever choice of the reference constant 0 is crucial. For instance, inserting the upper bound 0 = + into the estimate (8) with = 0 produces For fixed tolerance, the latter estimate implies that the necessary iteration count is proportional to the material contrast + ∕ − . Thus, no advantage is gained over the basic scheme. 19 In particular, due to the character of the optimum reference constant expressed as a geometric mean of the constants ± , a good guess for the lower bound − is critical. In contrast, a good upper bound + is usually inferred from the elastic stiffness of the material. To conclude this section, let us briefly comment on the case of porous microstructures, where the stress operator is identically zero on a non-trivial subdomain of the unit cell Q. Then, the stress operator is no longer strongly monotone, and the Giselsson-Boyd estimate (8) no longer applies. The general convergence statement discussed at the beginning of this section ensures convergence, at least in theory. However, this convergence behavior can be excruciatingly slow in practice, even leading to divergence in case of accumulated round-off errors. This is the reason for debates in the community concerning whether polarization schemes converge for porous materials or not. 23,24,39 Recently, the convergence behavior of primal algorithms for porous microstructures was analyzed. 49 It was found, both computationally and numerically, that the discretization scheme plays an important role in the stability of such numerical methods for porous microstructures. Furthermore, in the continuous (non-discretized) case, the convergence behavior of the basic scheme (and its descendants) is closely tied to a specific subspace of compatible strain fields, on which the stress operator is strongly monotone. The monotonicity constant, however, strongly depends on the geometry of the solid-pore interface, and is not readily accessible.
Without going into details, it is readily seen from the iterative scheme (5) that the strain field k+ 1 2 is always compatible, and even lives in the specific subspace on which the stress operator is strongly monotone. Unfortunately, in the ADMM scheme, the stress operator is only applied to the incompatible strain field e via = (e), and the non-degeneracy of the stress operator cannot be transferred to the compatible strains, invoking a rigorous convergence statement. However, it seems reasonable that, if the residual |L 2 is sufficiently small (and a stable discretization is used), then ADMM (5) will enter a regime of linear convergence. This statement can be made rigorous, see Liang et al. 53 We will further study this issue numerically in Section 4.
The plan for the remainder of this section is to discuss a number of strategies for adaptively estimating the reference constant k , fixing the damping factor k ≡ .

Finite material contrast
Let us start with the simplest case of (non-degenerate) linear elastic materials. 23,24,54 Denote by ± (x) the smallest and largest eigenvalue of the local stiffness tensor C(x) at the point x ∈ Q. Then, we may compute , For the case of non-linear materials at small strains, the following strategy proves useful. 28 Denote by k,± (x) the smallest and the largest eigenvalue of the material tangent evaluated for the current strain field e k at iteration k. Then, we may set and define k+1 = √ k,− k,+ . Typically, k,+ is independent of k, and k is a monotonically decreasing sequence, compare condition 3 in Section 2.1.
This approach comes with the advantage that, for inelastic materials featuring an elastic region, loadings at small strain lead to convergence rates in the polarization scheme like for the elastic case. Clearly, care has to be taken for materials whose material contrast goes to infinity for increased loading 28(sec.7) .

Residual balancing
A classical approach, introduced by He et al., 45 for adapting the step size in ADMM is residual balancing. For the case at hand (5) and unrelaxed ADMM = 1∕2, we define the primal and the dual residual The letters p and d were chosen to stand for primal and dual, respectively. In particular, no confusion with the stress polarization is intended. The quantity p k+1 (10) measures the feasibility condition ! = e in the primal form of the optimization problem (1), whereas the increment d k+1 (11) may be interpreted as quantifying feasibility for the Lagrangian dual of problem (1), see Boyd et al. 22(sec.3.3) In our case, please note the following. For a start, algebraic manipulations reveal the estimates and where Γ ≡ Γ k ∶ C k is independent of k . We refer to Appendix A for a derivation. Thus, the primal and the dual residual measure the compatibility and the strain boundary condition as well as the equilibrium and the stress boundary condition individually. Also, due to the identity the primal and dual residuals may be interpreted as a decomposition of the total (consistent) residual (22), to be discussed below in Section 3. He et al. 45 suggested that the primal and dual residuals should decay at the same rate in order to ensure fast convergence of ADMM, and to adapt the step size k if necessary. Some care has to be taken, as the primal residual (10) is a strain field and the dual residual (11) is a stress field. Due to the mismatch of dimensions, no direct comparison is possible. Rather, it is necessary to define a suitable interval of salient stress values [ min , max ] that is considered appropriate and to ensure the condition directly via an amplification factor A > 1 in terms of the adaptive rule The rationale behind this strategy is that, if the primal residual is large, increasing the penalty factor k should lead to a decrease of the term penalizing the difference − e. Typical values for parameters are A = 2, max = + and min = + ∕10 22(sec.3.4.1) . The resulting strategy is called residual balancing. Computing the primal and dual residuals may be integrated into Algorithm 1 without memory overhead.
The original residual-balancing strategy of He et al. 45 was developed for = 1∕2. We discuss a simple modification for k ≡ > 0. Then, defining the primal and dual residuals by the identity as well as the estimates (12) and (13) follow, see Appendix A. The residuals for > 0 may also be computed without additional memory overhead. Please note the explicit expression for the dual residual In particular, for ≡ 0, the identity d k+1 = k p k+1 holds, and the residual balancing strategy does not work.

Lorenz-Tran-Dinh updating
Lorenz and Tran-Dinh 46 studied adaptive step-size strategies for the Douglas-Rachford splitting and the ADMM. They provided a novel convergence analysis for operator-splitting methods of Douglas-Rachford type with adaptive parameter with absolutely summable increments. Based on a careful analysis of the linear case, they proposed a strategy which is valid for possibly multi-valued maximally monotone operators. For the iterative scheme (5) at hand, their suggestion takes the form of a secant stiffness. For linear elasticity, the proposed reference material k is guaranteed to lie within the minimum and the maximum eigenvalues of the local stiffnesses.

Barzilai-Borwein step size
Barzilai and Borwein 43 introduced an adaptive strategy for selecting the step size s k in gradient descent methods via the secant-type rule where angular brackets denote the inner product. The Barzilai-Borwein method turned out to be a powerful solution method for unconstrained smooth optimization, 55-57 despite its apparent simplicity. The Barzilai-Borwein (BB) step-size rule (17) was shown to be successful in FFT-based computational micromechanics, 39 as well.
For the ADMM, Xu et al. 47,48 studied step-size selection strategies of Barzilai-Borwein type. In this work, we consider a variation of the approach of Xu et al. More precisely, Xu et al. propose a BB strategy for the primal and the dual step size, and select the geometric mean of these two step sizes. For homogenization problems, it turns out to be more effective to rely upon a purely primal strategy. More precisely, the starting point of our discussion is the BB rule As the fields and e are related by the stress operator , we observe Thus, if the stress operator is − -monotone (6) and + -Lipschitz (7), the inequalities follow. Please note that Barzilai-Borwein 43 also introduced a second step-size rule, which, in our case, becomes The Cauchy-Schwarz inequality implies the ordering Thus, the second strategỹk +1 is more conservative than the first strategy (18). Quite often, the second strategy is slower in practice, and we will work with the first choice (18).
Please note that the bounds (19) are only effective for strongly convex stress operators. Indeed, if − = 0, for instance in the case of porous microstructures, the reference material may decay to zero, leading to an unstable, and eventually diverging scheme. For this purpose, we augment the rule (18) by a safeguard proposed by Xu et al. 47(sec.3.3) For ∈ (0, 1), we define This safeguard ensures that the iteration-dependent lower bound holds, which may be connected to the interface constant c I . 49(sec.2.2) Please note that the strategy (20) may be integrated into the ADMM Algorithm 1 without memory overhead. Last but not least let us comment on the selection of . We follow the suggestion of Xu et al. 47 and fix = 0.2.

The convergence criterion
In this section, we assume, for simplicity, that the reference materials C k commute with the projector P. Otherwise, the formulae become lengthy. 44 To assess convergence, we utilize the C k -and D k -weighted L 2 -norms Please notice that, in view of the identification of the ADMM (5) with a damped Eyre-Milton scheme, see Moulinec and Silva 24 , via the identification P k = k + C k ∶ e k , and as the identity k = (e k ) holds, it appears natural to consider the field e k as the primary quantity on which the ADMM scheme (5) iterates. In particular, we should quantify the deviation from the following equations Appendix B details the derivation of the identity Thus, we are led to the surprisingly simple result that the four different conditions (21) to be satisfied may be quantified by monitoring the quantity e k − k+ 1 2 only. Several remarks are in order.
1. For the basic scheme with mixed boundary conditions the identity holds. 28(eq.5.7) As each iterate k (k > 0) of the basic scheme (23) is kinematically compatible and satisfies the P-part of the strain boundary condition, the ADMM residual (22) appears as a natural extension of the residual (24) for the basic scheme. In fact, a convergence criterion implemented for the basic scheme may be reused for the ADMM, surprisingly. 2. For the pure ADMM case ( = 1∕2), Michel et al. 19(sec.4.1.2) suggested to monitor ||e k − k || and || k − ( k )||. If the stress operator is Lipschitz continuous with constant L, we observe Thus, monitoring the quantity k − ( k ) is not necessary and may be avoided for performance reasons. 3. For non-changing reference material C k ≡ C 0 , the "natural" convergence criterion for the damped Eyre-Milton iteration 28 with mixed boundary conditions, involving the operators where 0 ( ) ≡ C 0 ∶ , is motivated by the identity 28(eq.5.8) for the identification P k = k + C 0 ∶ e k with k = (e k ). Comparing to the expression for the ADMM (22), we notice the identity Thus, the residual for the ADMM is in fact identical to the residual for the damped Eyre-Milton method for non-changing reference material. Also, the established identity (22) may be considered as an extension of the identity (25), valid for polarization schemes, to non-stationary reference material.

Efficient implementation
The algorithm (5) for mixed boundary conditions is straightforward to implement on three strain-like fields, see Algorithm 1. Thus, the damping parameter , whose clever choice may significantly decrease the time to solution, does not lead to additional memory requirements compared to the original version of Michel et al. 19 The Algorithm 1 also indicates suitable instances when to update the reference material and when to compute the residual. Notice that no extra memory is required for evaluating the residual. For suitable initialization and stationary reference material, the Algorithm 1 produces identical iterates (and identical residuals) as Alg. 1 in Schneider et al. 28 However, the implementation of Algorithm 1 requires an additional field to be stored. Still, the original implementation 19 appears to be the most popular among practitioners of FFT-based methods, and is extended to Algorithm 1 with minimal effort. ← + C 0 ∶ ( − e) ⊳ Or as a byproduct in line 11 13: Update reference material C 0 14: until res < tol ⊳ Compare section∼3.1 15: return ( , ) To complete the section on the implementation, let us specify that we use a discretization on a staggered grid, 36 see Schneider 2(sec.2.5) for the efficient evaluation of the Γ-operator. Please note, however, that the implementation is suitable for alternative discretizations, as well, see the overview article. 2(sec.2)

Setup
The methods described in this article were integrated into an in-house FFT-based computational micromechanics code (written in Python with Cython extensions) which is parallelized using OpenMP. All examples in this section are solved up to a tolerance of tol = 10 −5 , measured in terms of the corresponding consistent convergence criteria. The run times refer to a desktop computer with six 3.70 GHz cores and 32 GB RAM. The different material parameters used in this article are summarized in Table 1, and we refer to Figure 1 for an impression on the considered microstructures.

A short-fiber reinforced composite
Our first example concerns a short-fiber reinforced polymer. Such a test is a standard benchmark in FFT-based computational micromechanics because it permits studying finite material contrast, but complex microstructure features. More precisely, we consider short glass fibers with a length of 200 μm and a diameter of 10 μm at 15% volume fraction inside a cubic volume element with edge length of 288 μm and an isotropic fiber orientation. The microstructure, see Figure 1A, was generated with the SAM algorithm 60 and comprises 229 fibers. We use a rasterized image with 192 3 voxels, that is, a mesh size h = 1.5 μm. We furnish the phases with the material parameters of the polymer matrix and the E-glass fibers listed in Table 1. As a warm-up, we performed simulations in the linear elastic setting, see Table 2, for 5% uni-axial extension in x-direction. We study the performance of ADMM with different relaxation factors , corresponding to the accelerated scheme of Michel et al. 19 for = 1∕2, the Monchiet-Bonnet recommendation 27 = 1∕4 and the Eyre-Milton scheme 25 for = 0. We observe that, for all adaptive choices of the reference material, the iteration count increases with , that is, in order to minimize the iteration count, setting = 0 is advised. In general, the optimum choice, see Section 2.2, which is available due to the finite material contrast, performs best. All other adaptive strategies are initialized with this theoretically optimum choice. Residual balancing, for instance, retains this optimum choice. Thus, its iteration counts are identical to those of the optimum choice (except for = 1∕2 where the reference material is changed once, and directly changed back after the next iteration). Please also keep in mind that residual balancing only works for ≠ 0.
In addition, we also considered the choice 0 = + ∕5 which performs well in Monchiet-Bonnet. 27( fig.4) . Indeed, the performance of the basic scheme and the ADMM will be quite similar if identical reference materials are chosen. Rather, the basic scheme is restricted by the condition 0 > + ∕2 to retain stability, 19 whereas no such limitation applies to the polarization schemes. Therefore, 0 = + ∕5 represents a rather aggressive choice for the reference material in the polarization scheme. For the linear elastic example at hand, this choice leads to slightly slower iteration counts than the Lorenz-Tran-Dinh scaling.
Both, the Lorenz-Tran-Dinh and the Barzilai-Borwein strategies perform worse than the optimum choice, but not dramatically so. Compared to the primal solvers, where we include the "golden standard", the linear CG, and the primal Barzilai-Borwein method, ADMM performs admirably. For optimum parameter choice, the required iteration count lies below CG. These findings are illustrated in Figure 2A which shows the residual versus current iteration for the different strategies and the damping factor = 0. We observe that the convergence rate of ADMM with optimum parameters matches the rate of linear CG. Both the primal Barzilai-Borwein method and the Barzilai-Borwein strategy for the ADMM, lead to a non-monotone behavior of the residual. However, these oscillations are less pronounced for the primal-dual version. The Lorenz-Tran-Dinh strategy is the fastest for low accuracy, up to 10 −3 , but experiences a decrease in convergence rate for higher accuracy. Next, we study the inelastic behavior. Using the material parameters of Table 1 for the von-Mises elastoplastic matrix with exponential-linear hardening, we subject the composite to 5% uni-axial strain loading, distributed over 50 load steps. The mean iteration counts and the run times are listed in Table 3. Except for residual balancing, increasing the damping factor also leads to an increase of the iteration count for the ADMM. The Lorenz-Tran-Dinh strategy is significantly faster than the theoretically optimum choice. Also, the Barzilai-Borwein strategy almost matches the speed of the optimum choice. Residual balancing, on the other hand, does not appear competitive. The constant strategy 0 = + ∕5 takes about twice as long as the optimum strategy.
Compared to the primal solvers, where we include the Barzilai-Borwein method, 39 nonlinear CG, 7 and Newton-CG 11,12 with the solver parameters identified in Wicht et al., 15 the ADMM compares favorably. Both Newton-CG (where we count, both, linear and nonlinear iterations) and nonlinear CG operate on a similar level as ADMM with the optimum choice. The (primal) Barzilai-Borwein method is significantly faster than the other primal methods, but cannot match the performance of ADMM with the Lorenz-Tran-Dinh scaling of the reference material. Indeed, the latter is twice as fast.
Let us take a look at the individual iteration counts per load step, see Figure 2B. Except for the very first steps, the Lorenz-Tran-Dinh scaling leads to the lowest iteration counts for every load step. This is consistent with the previous  observations, where this method was able to reduce the residual quickly. Combined with the affine extrapolation, this initially fast decrease of the residual leads to low iteration counts for this method. It is interesting to note that the "optimum choice" for the reference material in ADMM does not lead to the minimum iteration count, for instance exemplified by the peak at load step seven. The reasons are that, for a start, the optimum choice only realizes the optimum convergence rate. It is insensitive to the prefactor in front of the rate. In particular, combined to affine extrapolation, the optimum convergence rate is unable to come into play. Also, the optimum rate is only best if the bounds ± are actually sharp. In any case, studying adaptive choices for the reference material in ADMM can be extremely productive, also in the case where the supposedly optimum parameter choice is known.
The constant strategy 0 = + ∕5 performs well for the first ten load steps, but lags behind later on. Indeed, this strategy is fixed beforehand, and once and for all. In particular, the reference material cannot account for the increasing plastification for higher loading. All in all, this strategy (with = 0) is about 4.5 times slower than the Lorenz-Tran-Dinh scaling.

A closed-cell foam
In this section, we consider a closed-cell foam, see Figure 1B. Following Abdullahi et al., 61 the microstructure was generated by shrinking the cells of a Laguerre tessellation, where we use the Newton algorithm discussed in Kuhn et al. 62 to generate 27 Laguerre cells of equal volume. The resulting foam has a volume fraction of 11.6%, and was discretized by 128 3 voxels. Please note that we used a higher resolution for the visualization in Figure 1B. This foam example is rather challenging for FFT-based solvers. In fact, FFT-based methods are certainly not the optimum for treating microstructures with such a high degree of porosity (88.4%). Rather, the foam structure should be regarded as a benchmark showing the limitations of FFT-based methods, in general, and the ADMM, in particular. We furnish the solid foam with the isotropic and linear elastic material parameters of aluminum listed in Table 1. The iteration counts for 5% uni-axial extension in x-direction are listed in Table 4. We prescribed a maximum of 5000 iterations, and ADMM with residual balancing and = 1∕4 failed to converge within these limits. No damping ( = 0), leads to divergence of the solver.
The conjugate gradient method converges in about 400 iterations. The primal BB method requires more than 1000 iterations, and is closely matched by the primal-dual BB. For this example, the Lorenz-Tran-Dinh scaling requires significantly more iterations, roughly by a factor of 3, than the BB scaling. The 0 = + ∕5 strategy does not reach the prescribed tolerance within 5000 iterations.
Investigating the residuals more closely, see Figure 3C, we notice that the CG method does not converge monotonically, but rather shows a marked zig-zag pattern. This reflects the high condition number of the problem under consideration. Up to about 300 iterations, the CG residual is similar to BB's. After this point, CG significantly speeds up. All other considered ADMM schemes appear to have the same convergence rate, but with different prefactors. Please also note the severe oscillations in the residual of the primal BB scheme which cover three orders of magnitude.
We shall use the foam example to also investigate the influence of the discretization method and the appropriateness of the residual (22). Indeed, for finite material contrast, the various discretizations available for FFT-based computational micromechanics perform similarly well. Differences emerge for (highly) porous materials, like the open-cell foam. For this purpose, we re-ran the computational experiments with identical settings but different discretization schemes. We study the original discretization of Moulinec and Suquet, 3,4 which might be interpreted as an underintegrated Fourier-,Galerkin method, 63 and the discretization on a rotated staggered grid, 64,65 introduced to FFT-based micromechanics by Willot. 35 The residuals versus iteration are shown in Figure 3. We observe that none of the considered solution schemes, be it primal or primal-dual, converges within 5 000 iterations. For the Moulinec-Suquet discretization, all solvers fail to advance significantly below a residual of 10 −2 . More strikingly, the linear conjugate gradient method 5,6 stagnates at a residual above 10%. For Willot's discretization, the solvers fare slightly better, reaching residuals which are about an order of magnitude lower than the Moulinec-Suquet discretization. Also, a stagnating convergence behavior is observed for linear CG. This convergence behavior of the solution methods is in stark contrast to their staggered grid counterparts, see Figure 3C, which perform significantly better. Recall that such a discrepancy between the performance of primal solvers and different discretization schemes was already observed in Schneider et al. 36 To study the interplay of the discretization scheme and the residual (22), we monitor the effective axial stress per iteration in Figure 4. Strikingly, for the Moulinec-Suquet discretization, see Figures 4A, the effective stress decreases monotonically, and is expected to converge to zero as the iteration count goes to infinity. After 1 000 iterations, all the considered solvers give rise to different predicted effective stresses. Such an anomalous behavior is characteristic for the Moulinec-Suquet discretization applied to porous microstructures. 49 For Willot's discretization, see Figure 4B, the effective stresses appear to converge to a unique value, as is the case for the staggered grid, see Figure 4C.  Actually, for the ADMM solvers and Willot's discretization, the change in the effective properties is about an order of magnitude higher than for the staggered grid, reflecting the different magnitudes of the residuals, which we observed in Figure 3.
To sum up, for this highly porous example, the influence of the discretization scheme becomes clear. Only for a well-conditioned discretization, 49 a robust convergence behavior of FFT-based solvers can be expected. This holds for primal solvers as well as for primal-dual solvers. In this regard, the staggered grid discretization shines, as it appears to be the most robust for highly porous materials among the studied discretization schemes.
The computational results confirm that the "natural" residual (22) serves as an appropriate indicator for convergence (or lack thereof), independent of the considered discretization scheme. Yet, when used in conjunction with the staggered grid discretization, ADMM appears to reach its limitations. It converges, but cannot be considered competitive to the fastest primal solver.

A porous MMC
The last example concerns a Metal-Matrix Composite (MMC) with spherical pores (in blue), see Figure 1C. Apart from 29.5% spherical ceramic fillers, 0.5% spherical voids were included in the aluminum matrix. More precisely, 33 ceramic fillers of identical size were dispersed in the structure, together with 31 spherical pores with only 1∕4th of the radius of the ceramic fillers. The spheres were placed inside the structure by the mechanical contraction method. 66 The structure was discretized by 128 3 voxels, and we consider uni-axial extension up to 1% strain, distributed over 50 equidistant load steps. We consider the ceramic fillers to be linear elastic, and furnish the aluminum matrix with a von Mises elastoplastic constitutive model and a power-law hardening, see Segurado et al. 59 The specific model parameters are listed in Table 1. Due to the exponent in the power-law hardening, the smallest eigenvalue of the material tangent in the aluminum matrix goes to zero as the strain loading goes to infinity. Thus, this setup challenges the optimum choice 0 = √ − + even in the non-porous case, see Schneider et al. 28(sec.7.3) Of course, including spherical pores creates additional difficulties for the FFT-based solvers. The performance of the considered solvers is summarized in Table 5. We first take a look at the primal solvers. Newton-CG requires the highest iteration count. This is not unexpected, because the principal advantage of Newton-CG rests upon the fact that the iterations of the linear CG can be much faster than the nonlinear iterations. For the problem at hand, the constitutive models are not much more expensive than the linear elastic model. Therefore, Newton-CG and nonlinear CG lead to pretty much the same run time. In contrast, the primal Barzilai-Borwein scheme is much faster in terms of run time. This roots in the relative simplicity of the BB update. Indeed, it is little more than the basic scheme, and does not feature the more involved linear algebra operations necessary for nonlinear CG and Newton-CG. As the constitutive updates are so cheap, the linear algebra operations need to be accounted for.
For the ADMM, the Barzilai-Borwein scaling does not perform well, at all. Residual balancing with = 1∕4 leads to roughly the same iteration count as primal BB. However, the run time is much higher. Indeed, the linear algebra operations necessary for the ADMM significantly exceed those for primal BB. The Lorenz-Tran-Dinh scaling, on the other hand, slightly outperforms primal BB in terms of run time, and does so by requiring significantly less iterations. The iteration-independent choice 0 = + ∕5 is about 56% slower than the fastest solver.

CONCLUSIONS AND PERSPECTIVES
In this work, we studied strategies for choosing the reference material adaptively for the alternating direction method of multipliers introduced in the context of FFT-based computational micromechanics by Michel et al. 18,19 The conclusions are the following.
1. In the literature, it was not fully clear whether the accelerated scheme converges for porous materials or not. 23,24,39 The confusion may be attributed to different convergence criteria, a different choice of reference material, and usage of ill-conditioned discretization schemes.
In this work, we reported examples for a stable convergence process, and underlined these findings by theoretical arguments, see Section 2.1. More precisely, for a stable discretization scheme, for instance the staggered grid, 36 the damped ADMM Algorithm (5) converges for non-vanishing damping parameter . Under a geometric assumption on the solid-pore interface, 49 the convergence is linear.
Still, optimized primal solvers outperformed the considered primal-dual solvers. It remains to investigate whether an improved selection strategy for the reference material permits outperforming primal solvers for porous microstructures, as can be done for finite material contrast. 28,42 2. We devoted attention to adaptive choices of the reference material in the damped ADMM. Such adaptive choices were shown to represent memory-efficient alternatives to the Anderson-accelerated polarization methods, 42 sometimes even improving upon their performance. 3. For the ADMM, the residual-balancing strategy 45 is a classic, and widely used. We discussed the strategy for the problem at hand, see Section 2.3. In the process, we extended the strategy to variable damping, and identified the individual primal and dual residuals. It is interesting to note that, for the cell problem of micromechanics, the primal and dual residual vectors may be linearly combined to the consistent residual vector. This appears not to be the case in general. 22 Based on these results, we recommend this strategy. 5. We investigated a Barzilai-Borwein step size for the ADMM, 47 see Section 2.5. In contrast to the Barzilai-Borwein primal solver, 39 the BB strategy leads to significantly less oscillations in the associated residuals. Still, the strategy inherits the convergence speed from the primal setting, and is also recommended. 6. For primal FFT-based solvers, it was shown 7 that choosing the proper step size may be of secondary importance, as long as the momentum term is chosen in an adaptive manner. For the work at hand, we fixed the damping parameter . It might be of interest to investigate whether choosing the damping parameter adaptively can be beneficial, see Xu et al. 48 7. In this work, we considered convex and differentiable objective functions. The ADMM may be applied to non-differentiable objective functions, as well. 32,33 Furthermore, for computational homogenization at finite strains 10 , ADMM was applied successfully, 13,31 as well. Adaptive strategies for the reference material should be of interest for these scenarios. Similarly, treating non-convex problems like phase-field fracture on microstructures 9,67 by ADMM could be pursued.

ACKNOWLEDGMENTS
This work was supported by the German Research Foundation (DFG) within the International Research Training Group "Integrated engineering of continuous-discontinuous long fiber reinforced polymer structures" (GRK 2078). The anonymous reviewers deserve gratitude for providing constructive feedback.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A. IDENTITIES FOR PRIMAL AND DUAL RESIDUALS
In this section, we wish to show the estimates (12) and (13) for the primal and dual residuals (14), respectively, and the general algorithm (5) To derive the estimate for the primal residual (A1), we utilize the compatibility of k+ 1 2 k+ 1 compare Equation (B4), and write ) .
Taking norms on both sides and using the Helmholtz decomposition (B2) shows the estimate (A1). For the dual residual, the last line of the iteration (A4) implies In view of Equation (B4) in the form we may insert the latter identity into the first line (A4) to obtain that is, the identity holds. Taking norms on both sides and using the non-expansivity of Γ k ∶ C k yields the dual estimate (A2).

APPENDIX B. DETAILS ON THE CONVERGENCE CRITERION
We wish to show the identity Before coming to the details, let us record the following consequence of the Helmholtz decomposition in (linearized) elasticity for mixed boundary conditions 28(Apx.) valid for any square-integrable strain field, measuring (in order) equilibrium, compatibility and the magnitude of the mean, decomposed into the Pand Q-parts. Let us investigate the first line in the Algorithm (5) in more detail. Separating the homogeneous and non-homogeneous parts, we observe As Γ k ∶ C k is the C k -orthogonal projector onto kinematically compatible fields, it follows In order to assess the validity of the last three Equations 21, we need to monitor the quantity In view of the identity (B4), the latter quantity may be rewritten in the form To assess the equilibrium of the stress field k , comparing the Equations (B3) and (B4) yields the identity which might be rearranged into the form .
Inserting the latter identity into the formula (B5) yields Upon taking the (squared) C k -norm, the Helmholtz decomposition in the form (B2) implies the formula (B1).