Locally-synchronous, iterative solver for Fourier-based homogenization

We use the algebraic orthogonality of rotation-free and divergence-free fields in the Fourier space to derive the solution of a class of linear homogenization problems as the solution of a large linear system. The effective constitutive tensor constitutes only a small part of the solution vector. Therefore, we propose to use a synchronous and local iterative method that is capable to efficiently compute only a single component of the solution vector. If the convergence of the iterative solver is ensured, i.e., the system matrix is positive definite and diagonally dominant, it outperforms standard direct and iterative solvers that compute the complete solution. It has been found that for larger phase contrasts in the homogenization problem, the convergence is lost, and one needs to resort to other linear system solvers. Therefore, we discuss the linear system’s properties and the advantages as well as drawbacks of the presented homogenization approach.


Introduction
Homogenization subsumes methods and approaches that aim to provide the effective material properties of materials with microstructures. In the context of linear constitutive laws, the effective constitutive tensor L * relates the volume average of a gradient field g (strains, temperature gradient, charge density gradient, concentration gradient, etc.) to its thermodynamic conjugate divergence free field q (stresses, heat flux, electric current, concentration flow, etc.), Locally, the constitutive law is location dependent. Often, L * is approximated by L or a more elaborate mean field approach, some of which are summarized in Ref. [1]. More accurate results are obtained by solving the boundary value problem of a large virtual material sample-a representative volume element (RVE)-and extracting the average fields from the numerical solution. In the linear case, L * is then reconstructed by at most n such solutions, where n is the dimension of the vector space of g and q.
In the context of examining RVE realizations of microstructures, Fourier methods are an indispensable tool for various reasons. They are typically applied in different ways which are summarized in the following: • Microstructural data is often provided as pixel (2D) or voxel (3D) data, which can be efficiently transformed to its discrete Fourier representation. • The Fourier function basis satisfies the boundary conditions of the most common RVE problem of a periodically repeatable unit cell a priori. • The Fast Fourier Transform (FFT) allows to switch quickly between real and Fourier space, which allows for a fast evaluation of differential operators. • In their Fourier representation, the orthogonality of divergence-free and rotation-free fields (i.e., the Helmholtz decomposition of a vector field) can be exploited algebraically to derive a closed-form expression for the effective constitutive law.
That is to say, two distinct features-the availability of the fast transformation and the algebraic orthogonality of divergence-and rotation-free fields in Fourier space-are exploited differently, but both methods are summarized under the Fourier label. The first kind of Fourier (spectral) solver is a fixed point iteration scheme introduced by [2,3]. It solves a concrete RVE problem by rewriting the homogeneous partial differential equation (PDE) with non-constant coefficients that is to be solved. On the RVE domain the local material tensor L(x) = L 0 + L(x) is split and the L(x)-part is moved to the right-hand side, interpreting the new equation as an inhomogeneous PDE with constant coefficients. Solving this auxiliary problem successively corresponds to a fixed point iteration for the original problem. The easy transformation between real and Fourier spaces enables a fast evaluation of the right-hand side, i.e., an efficient iteration, by taking the derivatives in the Fourier space and applying the constitutive law involving L(x) in real space. Per iteration, one forward and backward transform are required. The effective properties are obtained as the linear relation between the homogeneous parts of the solution fields. For a complete identification of the effective properties, several RVE problems need to be solved (e.g., three average temperature gradients that are prescribed in orthogonal directions) or material symmetries need to be presumed. The main advantage of this method is that only function evaluations are needed, i.e., we do not explicitly need to solve a linear system. Therefore, this method is less memory-consuming compared to standard homogenization approaches, such that discretizations of 256 3 voxels are feasible [4]. Note that only the discretized fields of the microstructural data and the q, g, L-fields need to be kept in the computer's memory. Such schemes are usually referred to as spectral or FFT solvers. Variants of this approach can be found in the pertinent literature. Some effort went into reducing the already low memory requirements [5][6][7][8] and improving the iteration's convergence rate [9][10][11][12]. A more complete summary can be found in a recent paper by Wicht et al. [4]. In the remainder of this article, we refer to this method as "Moulinec-Suquet (MS) fixed point iteration".
The second kind of spectral solver exploits the Helmholtz decomposition, i.e., the orthogonality of gradient (rotationfree) and divergence-free fields, which is algebraic in Fourier space but differential in real space. The algebraic orthogonality in Fourier space allows to introduce projectors based on which a closed-form expression for the effective properties is derived. These methods rely on the linearity of the constitutive law, which is merged with the linear projection operations. To distinguish this method from a fixed-point iteration, we refer to it as the "projection method". An excellent account can be found in the textbook by G. Milton [13] in Sect. 12. This method is entirely formulated in the Fourier space. The fast Fourier transform is only needed once to set up the linear system in Fourier space. The latter is quite large, but needs to be solved only for a few variables. Approximate variants of this approach have been presented in Refs. [14][15][16][17][18]. At least to the authors' knowledge, a concise formulation of this approach has been first presented by Chen in Chap. 5 of Ref. [19]. She managed to solve such systems for discretizations d of the unit cell up to 11 3 voxels (data points). That is rather limited compared to today's standards, which is due to the memory intensive linear systems (see Sect. 5). It has complex coefficients, is unsymmetric, not sparse, and without a band structure. For the sake of clarity, the two approaches are compared in Table 1.
In the present contribution, the focus is on the second class of methods, which is by far the less popular of the two Fourier-based approaches. This is due to the much larger numerical effort, compared to the fixed-point method, because an unsymmetric non-sparse linear system with complex coefficients is solved for the tensorial auxiliary field (constitutive tensor, stiffness or conductivity), which has moreover more components than the usual gradient field (strains, temperature gradient) that is solved for in the fixedpoint iteration approach. However, we found that on the numerical side, some unused potential can be leveraged. Because only the homogeneous part of the auxiliary tensorial field is needed, and this is represented by the zeroth Fourier coefficient, only few solution components need to be computed. Thus, instead of applying a general purpose linear solver, a Jacobi iteration for the relevant solution components is used. Although the resulting scheme is similar to the fixed-point iteration method in the sense that the solution is iterated, it is still a different approach.
The article is divided into six sections including the introduction. Before we discussed the proposed methodology, the basic notation is explained in Sect. 2. The Fourier method is recalled in Sect. 3, while we describe the local iterative solver in Sect. 4. With this the theoretical foundation for our method is laid and we examine a specific example in Sect. 5. Here, also a comparison of the novel method and standard solvers in terms of their computational efficiency is provided. In Sect. 6, a brief summary of the most important findings is provided.

Notation
In this article, a direct notation is preferred. Vectors are denoted as bold minuscules, like p for the polarization, q for the heat flux or g for the temperature gradient. All tensorial quantities are represented w.r.t. the orthonormal basis e i when their components are considered. Locations in real space are denoted by x, for which we use Cartesian coor-dinates, i.e., x = x i e i . Second-order tensors are written as bold majuscules, such as L for the heat conduction tensor. The dyadic product and scalar contractions are denoted by (a ⊗ b ⊗ c) : (d ⊗ e) = (b · d)(c · e)a, with a dot being the usual scalar product between vectors. This notation is mostly needed for linear mappings between vectors, as in q = L · g. When the product is clear from the context, as in the last equation, the scalar dot is occasionally omitted. We further use the nabla operator ∇ and projection operators { }. The latter involves the projectors K = (k · k) −1 k ⊗ k and I − K, which project a vector onto its part parallel and perpendicular to k, where k is the wave vector of the Fourier base function.

Fourier method
In the current section, we present the fundamental ideas of the applied Fourier method. They can be found as well in the textbook of Milton [13] in Section 12. To this end, the Fourier series, the decomposition of periodic fields, the polarization problem, and the solution for effective material properties are discussed based on a generic problem and later for elastostatics and heat conduction as specific example problems.

Fourier series representation
For the sake of completeness, a brief introduction to Fourier series representations of functions is provided in the following. Let denote the scalar product of two functions, where the overline denotes the complex conjugate. It is easy to see that the functions form an orthonormal basis, i.e., where δ kl denotes the Kronecker(-delta) symbol and i= √ −1 is the imaginary unit. With respect to this basis, a sufficiently smooth function can be represented as an infinite series of f k (x) where the componentf k (also called amplitude or spectral coefficient) is given by the scalar product with f k (x) A generalization to functions defined in R 3 is obtained by specifying the direction, usually referred to as the plane wave vector k with integer components where k i are the components of k w.r.t. the orthonormal basis e i . The confinement to integer k i is due to the Fourier series expansion. In the next step, the plane wave vector k is decomposed into a unit direction vector n k and its magnitudek n k represents the direction of the plane wave, whilek is the wave number, which is connected to the frequency f and the wavelength λ, respectively, by The scalar product defined by Eq. (2) is then Again we see that where δ kg may be considered as a generalized Kronecker symbol With this, the approximation of a function defined in R 3 by a Fourier series becomes The components of the function in the Fourier space are obtained via a scalar product

Decomposition of periodic fields into three orthogonal parts
In elastostatics and Fourier heat conduction we have the following relations that govern the behaviour of a system, with the local balances, the constitutive laws and the integrability condition for the driving force of the process (listed from left to right): The same mathematical structure is found for Fick's diffusion, Ohmic resistance, or any process in which the flux of a conserved quantity is linear in the gradient of a field. Being concerned with a homogenization task, we removed the source terms from the balance equations, which are a heat source (or sink) in heat conduction problems and a volume specific force ρb and the inertia force density ρẍ in elastostatics. The reason for this is that these are not material properties that can be homogenized, but externally applied, material independent effects. For the sake of notational convenience, we drop the dependence on the location x of all fields. The differential constraints on ε and g are the so-called integrability con-ditions which ensure the existence of a field in case of the gradient g = ∇ f , where the interpretation of f depends on the physical problem and may be the temperature field, a concentration, or electric charge. In case of elasticity, the integrability conditions ensure the existence of a displacement field u for the strain field ε = (u ⊗ ∇ + ∇ ⊗ u)/2.
Next, we consider all periodic fields on a unit cell . Depending on the homogenization task, these are {q, L, g} or {σ , C, ε}. We proceed, for convenience, with the heat conduction problem. One can directly generalize the derivation given in the following to elastostatics. The divergence is obtained by the scalar product of f k with k Therefore, we can expand the Fourier series by the projector without affecting the divergence-part of the field With this definition we have f(x) · ∇ = f ·∇ (x) · ∇. An analogous calculation for the rotation gives This time, the part off k parallel to k can be removed due to the vector product of parallel vectors being zero. Thus, we can expand with I − K without affecting the rotational part of the field When normalizing k, one can see that the part for k = o is problematic. It represents the homogeneous mean value of the field, which is both divergence-and rotation-free. It needs to be addressed directly in the decomposition into three parts As can be seen, K and I − K filter the corresponding parts from the function without the imaginary unit. The imaginary unit appears upon taking the derivative, but we want to decompose the function itself. In real space, the orthogonality is easy to see as well with the permutation symbol ε i jk and Schwarz' theorem. The derivatives with respect to the coordinates x i are denoted by ,i . The Helmholtz-decomposition of f(x) into rotation-and divergence-free parts is not easily obtained in real space. But in the Fourier space, we just define three different Gammaoperators : These projectors extract the mean value ( 0 { }), the divergencefree ( ×∇ { }) and the rotation-free fluctuating ( ·∇ { }) parts of a field. Note that the application of the -operators is an algebraic operation, but not a simple matrix-product, and the application is not associative, i.e., A( B) = (A )B. Therefore, we denote the action of the operators by curly brackets as { }. It is clear that the three parts are orthogonal to each other w.r.t. the scalar product Elasticity The equations of elastostatics obey the same structure as the equations of heat conduction, cf. Eq. (15), which results in a similar decomposition of symmetric tensor fields into the homogeneous, divergence-and rotation-free parts: In the latter equation, one can easily identify the projection operators, where the projector that removes the diverging part of A involves the fourth order identity on symmetric second order tensors I. The strain tensor field ε(x) is rotation-free, hence its Fourier series can be written as while the stress field σ (x) is divergence-free and symmetric, hence its Fourier series can be written as

The polarization problem
We consider the same boundary value problem with an inhomogeneous conductivity L(x) and a homogeneous reference We apply homogeneous boundary conditions g on the entire surface. This makes the second problem trivial: all fields are homogeneous, moreover we have where g(x) denotes the fluctuation part of the g-field. Our interest is to find the linear mapping between the homogeneous parts which constitutes the implicit definition of the effective conductivity L * . We take the difference between Eqs. (31) and (32) which is more compactly written as Note that here only denotes a difference and should not be confused with the Laplace-operator. We move L 0 g(x) to the left-hand side where p(x) is defined as or in expanded form This is the so-called polarization problem. The name is most likely related to its first area of application which was electric polarization. Note that p(x) is neither divergence-nor rotation-free. This underlines its character as an auxiliary problem, which is not solved, but used to derive an explicit expression for L * . Taking the homogeneous part by applying 0 gives which can be brought into the form With an explicit expression for p, which is linear in g, one can identify the effective conductivity.

Solving for the effective conductivity
All fields, except the reference conductivity L 0 , the mean values q, p and g and the effective conductivity L * depend on the location x in real space or the wave vector k in Fourier space. For the sake of convenience, we drop these dependencies in the remainder of this section. Further, we use L 0 = l 0 I, which simplifies calculations. Apart from that, the choice of L 0 is not important for us, since we obtain a closed form expression for L * . For spectral solvers of the first type (see Sect. 1), the choice of L 0 is crucial for the convergence. Next, we solve symbolically for p with the aid of the -operators. We firstly apply ·∇ {L −1 0 { }} to Eq. (38): Since g = g + g ·∇ is rotation-free, ·∇ { } removes merely the homogeneous part of g. Since q is divergence-free and L 0 is a multiple of the identity tensor, we have ·∇ {L −1 0 q} = o. Thus, we obtain Abbreviating { } = ·∇ {L −1 0 { }}, this can be rearranged for g = g − {p} and inserted back into the Eq. (37) for p, For the sake of brevity, we define L − L 0 = L. This system of equations can be solved for p The inverse is not a matrix or tensor inverse, but the inverse action of the operator in the round brackets, which would be the integral with the operator's corresponding Green function. The inverse of the pure -operators clearly does not exist, since they project some part onto zero. But with the identity I inside the round bracket and the flexibility to choose l 0 as we please we can always circumvent singular operators. We now extract the mean value by applying 0 { } In the next step, we insert the result into Eq. (40) and factor out g as it does not depend on x Therefore, the effective conductivity is The action of the 0 -operator is to just take the volume average, which is usually denoted with . In the Fourier space, this is done by setting all Fourier coefficients to zero except the one that belongs to k = o. We can easily determine L * once we find the auxiliary field L , which has no direct physical interpretation. Instead of inverting the operator I+ L , we consider the solution of the linear system in the Fourier space: The second summand in the last equation represents a convolution. Because we have discrete periodic data in real space, the Fourier data is as well discrete and periodic. Hence, the summation is only carried out over one unit cell in the Fourier space, as depicted in Fig. 1.
A comparison of coefficients w.r.t. the function basis e ix·c in Eq. (53) requires b = c − k and a = c and yieldŝ which constitutes a linear system for the unknown Fourier coefficientsL c . One can see that this formulation has distinct mathematical advantages that can be directly exploited: The principal diagonal contains a sum ofL c and a product with coefficients that should decay at high frequencies.
Most important, we are only interested in the homogeneous part of the solution, i.e., theL o -component, see Eq. (49). This fact allows us to apply an efficient locally synchronous linear solver.
The presented system of equations is rarely set up and solved. Mostly, a series expansion for L * is used, see e.g., Chap. 14 in Ref. [13] or Sect. 11.2 in Ref. [20]. Then, the problem is addressed in the real space, which involves correlation functions of the microstructure and integrals over Green's functions. The latter integrals are divergent, hence a renormalization is necessary [21]. Even if this issue can be avoided, the nth order coefficients involve integrals over npoint correlation functions, which are quite tedious to obtain. Moreover, since one can control by the choice of L 0 whether a series alternates or converges from above or below, series expansions can give bounds for the effective properties.
As mentioned in the introduction (Sect. 1), other approaches solve for the gradient field g(x), which appears in an integral equation that has the same form as the Lippmann-Schwinger equation for quantum scattering. It is obtained by inserting Eq. (38) into Eq. (44). In this case, iterative methods to determine g(x) such as the alternating Real-Fourier space iterator [2,3] can be exploited. In this approach, several specific boundary problems must be solved if no symmetries are present. A variety of improvements and derivations of this method has been published, see e.g., Refs. [4,22].
In the next section, we will discuss a new type of iterative solver for systems of linear equations which will be used to determine L and ultimately L * = L 0 + L . This particular approach is capable of efficiently computing individual components of the entire solution vector which is in line with the proposed Fourier method, where only the homogeneous part of the solution is of interest to us.

Synchronous local algorithm for the iterative solution of linear systems of equations
In Sect. 3, it has been shown that analytical homogenization techniques lead to large linear systems despite the fact that only a few components of the result vector are of interest. Therefore, direct solvers are ill-suited to obtain solutions for such homogenization problems. In the following, we make use of a local solution algorithm that is capable of approximating only a few components of the solution vector. To achieve this goal, the sparsity of the coefficient (system) matrix has to be exploited. In the algorithm that is discussed in the following only the ith component of the solution to a system of linear equations is iteratively computed based on a Neumann series representation.

Motivation
The starting point for our discussion of the local iterative solver is a generic system of linear equations where A is a non-singular, real, square matrix of dimension n and b is the right-hand side vector of the system. Equation (55) can be solved by a wide variety of methods that can be classified into two major categories: (i) direct and (ii) indirect or iterative solvers. Typical representatives of the first class of methods are Gaussian elimination, LU-or Cholesky factorization to name just a few [23]. Iterative methods include stationary linear approaches (which will be discussed in more detail in Sect. 4.2) such as the Jacobi, Gauss-Seidel, or Richardson methods and gradient methods (conjugate gradients, preconditioned conjugate gradients, steepest descent methods) [24][25][26]. All of these methods have in common that the entire solution vector x is computed either exactly or in an approximate fashion.
To compute all components of vector x seems to be rather inefficient and a waste of computational resources if only a few components are actually of interest. It is easy to see that an iterative method that only computes selected components is much more efficient than a direct solution of the entire system of linear equations. For example, an LU-factorization of a fully populated matrix requires 2n 3 /3 floating-point operations (additions and multiplications), where n denotes the size of the system. The iteration procedure with fully populated matrices and residual requires 2n 2 per iteration [27]. We can expect, however, that much less than n iterations will be needed. Since only a few components of the solution vector are of interest in our analytical homogenization process, detailed in Sect. 3, we need to identify promising methods that are capable of computing individual components of the solution vector only. Therefore, the question regarding how this problem can be tackled. One possibility, based on simple stationary linear iterative methods, has been proposed by Lee et al. in a series of reports [27][28][29] and is discussed in the remainder of this section.

Stationary linear iterative methods
In order to keep the algorithm of the locally synchronous solution procedure as simple as possible, only stationary linear schemes are included in the discussion. In this contribution, we make use of different iterative solution schemes such as the (i) Jacobi, (ii) Gauss-Seidel, (iii) successive overrelaxation, and (iv) Richardson methods. The discussion of these methods closely follows Refs. [24][25][26]. The point of departure to derive algorithms for these methods is Eq. (55) which can be re-written as: where G denotes the iteration matrix, z the auxiliary solution vector, and the superscript (i) represents the ith step in the iteration. To ensure convergence of the results a necessary condition is i.e., the spectral radius of the iteration matrix must be less than unity 1 . This can be guaranteed if the system matrix A is positive-definite and/or diagonally dominant. Considering the matrix M, being introduced in Sect. 3, inequality (57) holds and consequently, iterative solvers are applicable. In general, each matrix A can be additively split into a sum of a strictly lower triangular matrix L, a strictly upper triangular matrix U, and a diagonal matrix D The solution to the system of equations (55) can be obtained iteratively by means of Eq. (56) which can be re-written using the matrixM which serves as a pre-conditioner of the system matrix A: with I being the unity matrix. Obviously, different choices for G and z are possible and will be addressed in the following, see Sects. 4.2.1 to 4.2.4.

Richardson method
The Richardson iteration technique is perceivably the simplest possible method and the prototype of any linear iteration. In this case, we choose the matrixM to be a multiple of the unity matrix, i.e., Therefore, the resulting iterative scheme is In this method, optimal convergence is achieved for a γ -value of where λ max and λ min denote the largest and smallest eigenvalues of the iteration matrix G R = I − 1/γ A. This choice minimizes the spectral radius ρ(G R ) of the iteration matrix and thus, maximizes the rate of convergence of the iterative method. If the coefficient matrix A is symmetric and positivedefinite and the Richardson method with an optimal choice 1 The spectral radius of a matrix B is defined as: ρ(B) = max{|λ| : λ ∈ σ (B)}, where σ (B) denotes the spectrum of B, i.e., the set of all eigenvalues of B.
of γ is used the spectral radius of the iteration matrix can be computed as indicating that ill-conditioned coefficient matrices (here κ( ) denotes the condition number) of linear systems of equations result in a slow convergence. This obviously limits the applicability for certain areas of application, e.g., large finite element systems.

Jacobi method
The Jacobi method is a simple iterative method that is used to compute the solution of linear systems of equations with a positive-definite and strictly diagonally dominant system matrix A, i.e., In this case, it is sufficient to choose the following set-up: By substituting Eq. (65) into Eq. (59), we obtain the iteration relation for the Jacobi method as: One important aspect to keep in mind is that the Jacobi method does not converge for every symmetric positivedefinite matrix. In order to improve the comparably poor convergence of the Jacobi method an additional weighting parameter ω can be introduced Considering the special case of a symmetric and positivedefinite coefficient matrix A, it can be shown that the weighted Jacobi method converges if [24] and the optimal choice of ω is provided for whereλ max andλ min are the largest and smallest eigenvalues of the matrix D −1 A, respectively. A typical choice for the weighting parameter ω is 3/2 which safes effort in computing the eigenvalues of the system. At this point we can derive the iteration expression as

Gauss-Seidel method
The idea of the Gauss-Seidel iterative solver is closely related to that of the Jacobi method (see Sect. 4.2.2). In contrast to the Jacobi method, each component of the new iterate depends on all previously computed components. Therefore, the following expression is used to set up the iteration: As in the previous section, Eq. (71) is substituted into Eq. (59) and thus the solution for each iteration can be computed as Note that the Gauss-Seidel method is applicable to linear systems if the coefficient matrix A is either strictly diagonally dominant or symmetric positive-definite. This method generally converges faster than the conventional Jacobi method, i.e., without weighting procedure, as presented in Sect. 4.2.2.

Successive over-relaxation
The successive over-relaxation (SOR) method can be directly derived from the Gauss-Seidel scheme, see Sect. 4.2.3, by extrapolation. To this end, the weighted average of the solution of the previous iteration step x (i) SOR and the solution of the Gauss-Seidel scheme for the current iteration x is taken. To derive the iteration sequence for the SOR scheme, a different splitting of the coefficient matrix A compared to the one introduced in Sect. 4.2 is used. Instead, overrelaxation is based on In the next step, Eq. (74) is substituted into Eq. (55) and we obtain the following relation to compute the solution iteratively We can reach the same solution by definingM as At this point, the question regarding the optimal choice of the extrapolation/relaxation factor ω remains. It can be shown that the optimal value is ω SOR,opt = 2 where μ is the spectral radius of the iteration matrix of the Jacobi method G J μ = ρ(G J ) and It is straightforward to prove that for ω = 1 the SOR is identical to the Gauss-Seidel method. Depending on the problem that is solved, the SOR method might be significantly faster than the Gauss-Seidel method. This has, however, not been observed for the analytical homogenization which is our main concern. An extension of the SOR method is the symmetric successive over-relaxation (SSOR) scheme requiring a symmetric coefficient matrix A. One step in this approach consists of a standard (forward) SOR computation followed by a backward one. To this end, the approach to splitting the coefficient matrix A presented in Eq. (74) is adjusted and therefore,M is now defined as Thus, the backward SOR step takes the following form The equations for the SSOR method can now be derived by substituting the solution obtained with Eq. (75) into the righthand side of Eq. (81). The final iteration expression can be set up using The iteration relation for the SSOR scheme is only given in the form of Eq. (59) due to the length of the resulting expression A suitable value for ω is still computed by means of Eq. (77). Nowadays, the SSOR scheme is primarily used as a preconditioner for other iterative schemes that are suitable for symmetric matrices as its convergence is usually slower than that of the standard SOR approach. Note that the optimal choice of the iterative method, or in other words, the optimal choice of the iteration matrix G and the auxiliary solution vector z are highly depend on the problem under consideration. Thus, for each method different properties of the system matrix A can be expected and these need to be accounted for.

Locally synchronous approximation of a single component of the solution vector
The synchronous algorithm is a modification of standard stationary linear iterative methods using an observation which allows the vectors involved in the computation to remain sparse when the iteration matrix G is sparse (with a maximum of d nonzero components; d n), 1/(1 − ||G|| 2 ) is sufficiently small, and n (number of equations) is large [27]. As solving a system of linear equations is fundamentally a problem involving the full matrix, it is not a trivial task to design an approach that extracts only a single component. The idea used in this contribution for the analytical (Fourier-based) homogenization is based on the works of Lee et al. [27][28][29]. Without loss of generality (considering the intended applications), the following discussions are limited to (i) positive definite and/or (ii) strictly diagonally dominant coefficient matrices A. In the first case, it is ensured that holds, while in the latter case is valid for symmetric iteration matrices.

Neumann series representation
All the iterative methods, introduced in Sects. 4.2.1 to 4.2.4, have in common that their primary goal is to approximate the leading term of the Neumann series As mentioned before convergence of this series is only ensured if the spectral radius of the iteration matrix ρ(G) is below one. The smaller this value is the faster the iteration converges, i.e., less steps are needed to fall below a prescribed error threshold. Different possible choices of G and z have been discussed above and can now be used to construct a local algorithm that does not solve for the entire vector x but only for a few components x i that are of interest due to various reasons depending on the intended area of application.

Extraction of a single component
The where the subscript (k) denotes the index of the component.
In the next step, we left-multiply Eq. (86) by e T i and obtain Thus, the vector z is only multiplied once at the end of the summation process. Since the Neumann series is an infinite series it has to be truncated after a certain number of steps t and therefore, the equation is split into two parts (89)

Algorithm
The complete algorithm including the initialization stage and the update procedure are summarized in Algorithm 1. For the sake of a compact notation, we introduce the residual vector r (t) and the estimate vector Thus, the iteration expression is Algorithm 1 computes the estimate according tô and thus, the error is given by That is also the reason why r (t) is coined residual vector as indicated above. Generally, it is observed that the method terminates within at most steps (iterations), and achieves an estimate such that holds, where is a prescribed error threshold. For us, the key property to exploit is the sparsity of the linear system, in conjunction with the damping effect of the iteration matrix on the residual: The initial iteration vector is zero everywhere except for the desired solution component (see Algorithm 1), for which it is one. Hence, a large number of products does not need an explicit evaluation. However, during the iterative process, the vector gets populated, hence more products need to be carried out. Depending on the Algorithm 1 Algorithm for the locally synchronous iterative solver. while ||r (t) || 2 ≥ do 9: p (t+1) = p (t) + r (t) ; 10: r (t+1) = G T r (t) ; 11: t = t + 1; 12: end while 13: returnx i = p (t)T z; x i : Approximation of the ith component of the solution vector 14: end function iteration matrix, the new entries are damped, or propagate through the residual vector. It turns out that the investigated linear system, resulting from the analytical homogenization approach discussed in Sect. 3, is appropriate for this iterative solver. When the iteration matrix G is sufficiently sparse and its spectral properties are well-behaved, i.e., d = O(1) and −1/ln(||G|| 2 ) = O(1), then we achieve an approximation for x i in constant time with respect to the size of the coefficient matrix for large n. For a more detailed discussion of the properties of this algorithm the reader is referred to the original work by Lee at al. [27][28][29].
The Richardson and Jacobi methods are particularly suitable for the local iterative solver if the coefficient matrix A is sparse as the iteration matrix G inherits this property. This will also lead to sparse vector z. Note that the approximation of the leading terms of the Neumann series will be at least as dense as z resulting in an increasing computational burden if the vector is fully populated or fill-in takes place. According to Ref. [27] bounds on the convergence rates and convergence rate guarantees can be provided if ||G|| 2 < 1 holds.
In order to apply the different iterative algorithms in conjunction with the concept of the locally synchronous solver, we only need to exchange the definitions of the auxiliary iteration variablesM and z in the procedure given in Algorithm 1. For the sake of convenience, the definitions of these variables are listed in Table 2 including a reference to the corresponding sections and equations.

Results
In this section we assess the performance of the presented locally synchronous iterative solver (see Sect. 4.3) in comparison to the default direct and Krylov subspace iteration methods provided in the computer algebra system (CAS) Mathematica.

Heat conduction problem
Since micro-structural data are generally provided on a Cartesian lattice, the linear system for L (see Eq. 54) is easily obtained with the help of the discrete Fourier transform (DFT). It provides as many data points in the Fourier space as are used in real space. As an example we consider a twophase material. We seek the effective heat conductivity L * that relates the effective heat flux q to the temperature gradient g: q = −L * · g. (97) Locally, we have in phases 1 and 2, respectively. The second law of thermodynamics requires the heat to flow from regions with higher temperature to those with lower temperature, hence the inequality q(x) · g(x) < 0 holds everywhere. We incorporate this negativity in the above equations explicitly, such that we can work with positive definite conductivities L i and thus, a positive definite L * .

Verification problem: Laminate structure
To check the implementation, a simple reference problem, i.e., a laminate structure with the isotropic local conductivities has been considered. This example has been chosen as a closed-form solution is available. Parallel to the plane the laminate the effective conductivity is the arithmetic mean value (Voigt average) of the conductivities l 1 and l 2 , while in other planes perpendicular to the lamination plane the effective conductivity is the harmonic mean (Reuss average) of l 1 and l 2 where the lamination direction is e 3 . The volume fractions v 1 + v 2 = 1 correspond to the plies' thicknesses.
As reference conductivity we choose l 0 = 1.5. The solution quality only depends on the ability of the discretization to capture the actual volume fractions. In cases where discretization exactly resolve the material distribution, e.g.,

Complex microstructure
As the second example, we consider a microstructure whose effective conductivity has to be computed numerically as no analytical solution can be derived. To examine the numerical properties of the worst possible case, we remove any periodicity, material symmetry, coordinate independence (as in the laminate) or symmetric alignment. The phase distribution is generated by periodically repeating an ellipsoid that is larger than the unit cell, see Fig. 4. The ellipsoid which is cut at the boundaries of the unit cell (dimensions: −1/2 ≤ x ≤ 1/2, 1/2 ≤ y ≤ 1/2, 1/2 ≤ z ≤ 1/2.) is defined by the following expression The ellipsoid penetrates the neighbouring cells, where it overlaps with copies of itself. In this way, we obtain a complex microstructure that has two interpenetrating phases.
No rotational symmetry is found in the arrangement of the microstructure, but still a cubic translational symmetry exists, which is unavoidable due to the periodicity of the Fourier base functions. The base vectors e i are parallel to the edges of the the cubic unit cell. A discretized version of this microstructure is depicted along with the Fourier coefficients in Fig. 5. The local conductivities are both orthotropic and not aligned with the cubic periodicity frame of the unit cell. This setup is the most general case that can be considered in the given framework: For second order tensors, the smallest material symmetry group is the orthotropic group, and the microstructure has no rotational symmetry.
The reference heat conduction is chosen as L 0 = l 0 I with l 0 = 3. This value is relatively close to the entries of L 1 and L 2 . If one uses series expansions to estimate L * , the choice of L 0 controls the rate of convergence. Here, the choice of l 0  LinearSolve[] 2 with the default behaviour, namely a direct solver, leads to significant computational costs. The solver times for different d-values ranging from 5 to 11 are compared in Fig. 6, while the corresponding convergence of the solution is visualized in Fig. 7.
From the first diagram we can infer that a solution time of approximately 20 minutes is needed for d = 11 when employing the internal direct solver, while only 8 minutes 2) requires approximately 20 seconds for this problem which is more than 20 times faster. For d > 11, the internal solvers run into memory issues, which happens only at d > 14 for the locally synchronous solver. The system matrix for a discretization of d generally exhibits a sparsity of approximately 45%, i.e., 45% of the matrix components are nonzero, and therefore has 0.45(6d 3 ) 2 entries. Even using sparse data formats, a fullfledged CAS is not an appropriate environment and therefore, tools that are more specialized on numeric computations are required. Nevertheless, the potential of the local solver becomes apparent.

Properties of the linear system
A typical system matrix plot is given in Fig. 8. The matrix is unsymmetric and has approximately 40% of nonzero entries, i.e., it is far from being as sparse as the usual finite element or finite difference stiffness matrices. We found the sparsity to depend on the microstructure, the phases' properties, the choice of the reference conductivity l 0 and the discretization. Two extreme cases for the laminate microstructure are illustrated in Figs. 9 and 10, respectively, where the sparsity depends on the parity of d. Note that the right-hand side vector has approximately 94.4% of nonzero entries.

Comparison to other linear systems
We have applied the same approach for finite element problems, where we tried to solve only for a specific degree of freedom (DOF). It turned out that none of the methods mentioned in Sect. 4 yielded a versatile iterative method that is comparable in terms of performance to standard direct solver. In our studies we also included different pre-conditioning techniques (diagonal pre-conditioner; column pre-conditioner [30]; upper co-diagonal pre-conditioner [31], maximum upper triangular matrix pre-conditioner [32]; incomplete Cholesky decomposition) and algorithms to make the coefficient matrices (strictly) diagonally dominant [33] to no avail. The top priority for all algorithms was placed on a simple and efficient implementation. While the incomplete Cholesky pre-conditioner is comparably effective, it is too costly to be useful in the proposed method. Additionally, minimizing the bandwidth of the system matrix using the reverse Cuthill-McKee method does not signif-  icantly improve the overall performance. We believe that the cause for this is that the solutions of boundary value problems require the boundary DOFs to communicate, in a sense. Boundary conditions on one side of the computational domain affect the DOFs on the other side, hence the information needs to propagate through the linear system. Unless the preconditioning is a quasi-inversion of the linear system, the locally synchronous iterative method cannot easily cope with such problems. Here, more sophisticated iterative solvers such as the biconjugate gradient stabilized method BiCGSTAB or multi-grid solver have to be employed and coupled with the local iterations (if possible).
On the other hand, in our application to homogenization based on Fourier methods we solve a linear system for Fourier coefficients, but are only interested in the leading order term. therefore, the higher the frequency gets, the smaller is the actual contribution of that coefficient. Hence, we observe a rapid damping of the contribution and the convergence of the local iterative method is ensured.

Limitations
At this point, important imitations of the presented approach need to mentioned. This facilitates a general understanding of the locally synchronous method and enables one to choose suitable areas of application. There are basically two facts that need to be considered: 1. Analogous to the MS fixed point solver, high phase contrasts pose a problem for the convergence. We managed to obtain converged solutions for phase contrasts up to 10 by adjusting l 0 . When the phase contrasts becomes larger, the optimal l 0 leaves some eigenvalues close to 1 in the spectrum of the iteration matrix, see Fig. 11, which prevents a convergence of the iterative solver. This is also shown in Fig. 12 for the lamellar structure. For a phase contrast of l 2 /l 1 = 8 and l 0 = l 1 , the largest eigenvalue is already larger than 1. In the special case of the laminate, the iteration still converges because no residuum appears in direction of the largest eigenvalue during the iteration. However, in a more general setting, convergence will not be observed.
Adjusting the weight ω of the Jacobi-preconditioning can accelerate the convergence, but does not lead to convergence when the standard Jacobi method with ω = 1 does not converge. It should be noted that this really is a problem of the solver: The linear system solution can be obtained for high phase contrasts with another linear solver. This flexibility is not available in case of the first kind of spectral solver. 2. The proposed approach is very memory intensive, which needs to be reduced in order to become a viable alternative to established numerical methods. It appears that the system matrix has a very specific structure, which may allow for an efficient compression. However, it is still questionable whether this is worth the effort.

Summary
Exploiting the orthogonality of rotation-and divergence-free fields algebraically in the Fourier space, we reformulate the difference between the effective conductivity and the reference conductivity as the solution of a large linear system in Fourier space, of which only the homogeneous solution components are needed, i.e., the zeroth Fourier coefficients. Therefore, instead of using a general purpose linear solver we employ a Jacobi iteration only for the solution required solution components. This enables us to apply standard methods to improve the convergence of linear iterative solvers. It turns out that for small phase contrasts, the method works very well and outperforms solvers that give the complete solution vector by several orders of magnitude. However, the iteration matrix has only an advantageous eigenvalue spectrum (which is direly needed for a fast convergence) for low phase contrasts and when an appropriate reference medium is chosen. For phase contrasts above 10, the solution is attainable by other means: either a more sophisticated pre-conditioning is needed or one has to resort to standard solvers. Unfortunately, while the solver is very fast, the linear system has huge memory requirements, such that its construction is not feasible for useful discretizations. The system matrix is unsymmetric, has complex coefficients, no band structure, and is with approximately 45% of nonzero entries not sparse. This leads to roughly 0.45 × (k × d 3 ) 2 complex matrix coefficients, where k is the number of components in the constitutive tensor, e.g., k = 6 for heat conduction problems, or k = 21 in linear elasticity. One can see that the memory requirements are excessive even for small d. Another limitation of the projector approach is that linear projections and linear constitutive laws are combined, such that an extension to nonlinear problems appears impossible. The above approach differs somewhat from the more usual spectral solvers introduced by Moulinec and Suquet [2], which define a fixed-point iteration for solving specific boundary value problems, from which then the volume average and the effective constitutive tensor is extracted. They solve for a lower-dimensional field, e.g., the strains or temperature gradient instead for the constitutive tensor, and without the need to set up the linear system, which is why this scheme is very memory efficient. The method's good performance relies on the fast Fourier transform that is needed in every iteration step. By linearizing constitutive laws, it is possible to consider nonlinear processes incrementally.
Our approach cannot compete with this method, but is rather complementary. For example, in the Moulinec and Suquet approach, only the volume averages of the solution fields are needed for identifying the effective constitutive tensor. It may be possible to combine the approach to solve only for the volume average in both methods.
Independently, the possibility to solve only for specific variables in linear systems that appear in engineering applications is interesting all by itself. We have shown that some performance gain may be achieved for a particular engineering problem in comparison to approaching the entire linear system directly. Our approach cannot compete with established methods for this particular problem, but it may be used in other engineering contexts. An application could be to calculate results in finite element calculations only at specific, critical locations. One example is the acoustic analysis of noise generation, e.g., engine sound, where the resulting sound pressure level is commonly only needed at few locations in the vicinity of the radiator [34]. Currently, linear solvers for specific variables seem to be employed in problems from information science only [27,35].