Analysis of a Fast Fourier Transform Based Method for Modeling of Heterogeneous Materials
Abstract
The focus of this paper is on the analysis of the Conjugate Gradient method applied to a nonsymmetric system of linear equations, arising from a Fast Fourier Transformbased homogenization method due to Moulinec and Suquet [1]. Convergence of the method is proven by exploiting a certain projection operator reflecting physics of the underlying problem. These results are supported by a numerical example, demonstrating significant improvement of the Conjugate Gradientbased scheme over the original MoulinecSuquet algorithm.
Keywords:
Homogenization, Fast Fourier Transform, Conjugate Gradients1 Introduction
The last decade has witnessed a rapid development in advanced experimental techniques and modeling tools for microstructural characterization, typically provided in the form of pixel or voxelbased geometry. Such data now allow for the design of bottomup predictive models of the overall behavior for a wide range of engineering materials. Of course, such step necessitates the development of specialized algorithms, capable of handling largescale voxelbased data in an efficient manner. In the engineering community, perhaps the most successful solver meeting these criteria was proposed by Moulinec and Suquet in [1]. The algorithm is based on the Neumann series expansion of the inverse of an operator arising in the associated LippmannSchwinger equation and exploits the Fast Fourier Transform to evaluate the action of the operator efficiently for voxelbased data. In our recent work [2], we have offered a new approach to the MoulinecSuquet scheme, by exploiting the trigonometric collocation method due to Saranen and Vainikko [3]. Here, the LippmanSchwinger equation is projected to a space of trigonometric polynomials to yield a nonsymmetric system of linear equations, see Section 2 below. Quite surprisingly, numerical experiments revealed that the system can be efficiently solved using the standard Conjugate Gradient algorithm. The analysis of this phenomenon, as presented in Section 3, is at the heart of this contribution. The obtained results are further supported by a numerical example in Section 4 and summarized in Section 5.
The following notation is used throughout the paper. Symbols , and denote scalar, vector and secondorder tensor quantities, respectively, with Greek subscripts used when referring to the corresponding components, e.g. . The outer product of two vectors is denoted as , whereas or represents the single contraction between vectors (or tensors). A multiindex notation is employed, in which with represents and abbreviates . Block matrices are denoted by capital letters typeset in a bold serif font, e.g. , and the superscript and subscript indexes are used to refer to the components, such that with
Submatrices of are denoted as
for and . Analogously, the block vectors are denoted by lower case letters, e.g. and the matrixbyvector multiplication is defined as
(1) 
with and .
2 Problem setting
Consider a composite material represented by a periodic unit cell
In the context of linear electrostatics, the associated unit cell problem reads as
(2) 
where is a periodic vectorial electric field, denotes the corresponding vector of electric current and is a secondorder positivedefinite tensor of electric conductivity. In addition, the field is subject to a constraint
(3) 
where denotes the dimensional measure of and a prescribed macroscopic electric field.
The original problem (2)–(3) is then equivalent to the periodic LippmannSchwinger integral equation, formally written as
(4) 
where denotes a homogeneous reference medium. The operator is derived from the Green’s function of the problem (2)–(3) with and can be simply expressed in the Fourier space
(5) 
Operator stands for the Fourier coefficient of for the th frequency given by
(6) 
”” is the imaginary unit (). We refer to [2, 4] for additional details. Note that the linear electrostatics serves here as a model problem; the framework can be directly extended to e.g. elasticity [5], (visco)plasticity [6] or to multiferroics [7].
2.1 Discretization via trigonometric collocation
The numerical solution of the LippmannSchwinger equation is based on a discretization of a unit cell into a regular periodic grid with nodal points and grid spacings . The searched field in (4) is approximated by a trigonometric polynomial in the form (cf. [3, Chapter 10])
(7) 
where designates the Fourier coefficients defined in (6). Notice that the trigonometrical polynomials are uniquely determined by a regular grid data, which makes them wellsuited to problems with pixel or voxelbased computations.
The trigonometric collocation method is based on the projection of the LippmannSchwinger equation (4) onto the space of the trigonometric polynomials
(8) 
leading to a to linear system in the form, cf. [2]
(9) 
where is the unknown vector, is the identity matrix, expressed as the product of the Kronecker delta functions and , and .
All the matrices in (9) exhibit a blockdiagonal structure. In particular,
with , and . The matrix implements the Discrete Fourier Transform and is defined as
(10) 
with representing the inverse transform.
It follows from Eq. (1) that the cost of multiplication by is dominated by the action of and , which can be performed in operations by the Fast Fourier Transform techniques. This makes the system (9) wellsuited for applying some iterative solution technique. In particular, the original Fast Fourier Transformbased Homogenization scheme formulated by Moulinec and Suquet in [1] is based on the Neumann expansion of the matrix inverse , so as to yield the th iterate in the form
(11) 
As indicated earlier, our numerical experiments [2] suggest that the system can be efficiently solved using the Conjugate Gradient method, despite the nonsymmetry of evident from (9). This observation is studied in more detail in the next Section.
3 Solution by the Conjugate Gradient method
We start our analysis with recasting the system (9) into a more convenient form, by employing a certain operator and the associated subspace introduced later. Note that for simplicity, the reference conductivity is taken as .
Definition 1
Given , we define operator and associated subspace as
Lemma 1
The operator is an orthogonal projection.
Proof
First, we will prove that is projection, i.e. . Since is a unitary matrix, it is easy to see that
(12) 
Hence, in view of the blockdiagonal character of , it it sufficient to prove the projection property of submatrices only. This follows using a simple algebra, recall Eq. (5):
The orthogonality of now follows from
according to a wellknown result of linear algebra, e.g. Proposition 1.8 in [8]. ∎
Remark 1
It follows from the previous results that the subspace collects the nonzero coefficients of trigonometric polynomials with zero rotation, which represent admissible solutions to the unit cell problem defined by (2). Note that the orthogonal space contains the trigonometric representation of constant fields, cf. [4, Section 12.7].
Lemma 2
The solution to the linear system (9) admits the decomposition , with satisfying
(13) 
Proof
With these auxiliary results in hand, we are in the position to present our main result.
Proposition 1
The nonsymmetric system of linear equations (9) is solvable by the Conjugate Gradient method for an initial vector with . Moreover, the sequence of iterates is independent of the parameter .
Proof (outline)
It follows from Lemma 2 that the solution to (9) admits yet another, optimizationbased, characterization in the form
(15) 
The residual corresponding to the initial vector equals to
It can be verified that the subspace is invariant, thus . Therefore, the Krylov subspace
for arbitrary . This implies that the residual and the Conjugate Gradient search direction at the th iteration satisfy and . Since is symmetric and positivedefinite on , the convergence of algorithm now follows from standard arguments, e.g. Theorem 6.6 in [8]. Observe that different choices of generate identical Krylov subspaces, thus the sequence of iterates is independent of .∎
Remark 2
Note that it is possible to show, using direct calculations based on the projection properties of , that the Biconjugate Gradient algorithm produces exactly the same sequence of vectors as the Conjugate Gradient method, see [9].
4 Numerical example
To support our theoretical results, we consider a threedimensional model problem of electric conduction in a cubic periodic unit cell , representing a twophase medium with spherical inclusions of volume fraction. The conductivity parameters are defined as
where denotes the contrast of phase conductivities. We consider the macroscopic field and discretize the unit cell with nodes^{1}^{1}1In particular, was taken consequently as and leading up to unknowns. The conductivity of the homogeneous reference medium is parametrized as
(16) 
where delivers the optimal convergence of the original MoulinecSuquet FastFourier Transformbased Homogenization () algorithm [1].
We first investigate the sensitivity of Conjugate Gradient () algorithm to the choice of reference medium. The results appear in Fig. 1, plotting the relative number of iterations for against the conductivity of the reference medium parametrized by , recall Eq. (16). As expected, solver achieve a significant improvement over method as it requires about iterations of for a mildlycontrasted composite down to for . The minor differences visible especially for can be therefore attributed to accumulation of roundoff errors. These observations fully confirm our theoretical results presented earlier in Section 3.
In Fig. 1, we present the total computational time^{2}^{2}2The problem was solved with a Matlab^{®} inhouse code on a machine Intel^{®} Core™2 Duo 3 GHz CPU, 3.28 GB computing memory with Debian linux 5.0 operating system. as a function of the number of degrees of freedom and the phase ratio . The results confirm that the computational times scales linearly with the increasing number of degrees of freedom for both schemes for a fixed [2]. The ratio of the computational time for CG and FFTH algorithms remains almost constant, which indicates that the cost of a single iteration of and method is comparable.
In addition, the memory requirements of both schemes are also comparable. This aspect represents the major advantage of the shortrecurrence CGbased scheme over alternative schemes for nonsymmetric systems, such as GMRES. Finally note that finer discretizations can be treated by a straightforward parallel implementation.
5 Conclusions
In this work, we have proven the convergence of Conjugate Gradient method for a nonsymmetric system of linear equations arising from periodic unit cell homogenization problem and confirmed it by numerical experiment. The important conclusions to be pointed out are as follows:

The success of the Conjugate Gradient method follows from the projection properties of operator introduced in Definition 1, which reflect the structure of the underlying physical problem.

Contrary to all available extensions of the scheme, the performance of the Conjugate Gradientbased method is independent of the choice of reference medium. This offers an important starting point for further improvements of the method.
Apart from the already mentioned parallelization, performance of the scheme can further be improved by a suitable preconditioning procedure. This topic is currently under investigation.
Acknowledgments
This work was supported by the Czech Science Foundation, through projects No. GAČR 103/09/1748, No. GAČR 103/09/P490, No. GAČR 201/09/1544, and by the Grant Agency of the Czech Technical University in Prague through project No. SGS10/124/OHK1/2T/11.
References
 [1] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, Comptes rendus de l’Académie des sciences. Série II, Mécanique, physique, chimie, astronomie 318 (11) (1994) 1417–1423.
 [2] J. Zeman, J. Vondřejc, J. Novák, I. Marek, Accelerating a FFTbased solver for numerical homogenization of periodic media by conjugate gradients, Journal of Computational Physics 229 (21) (2010) 8065–8071. arXiv:1004.1122.
 [3] J. Saranen, G. Vainikko, Periodic Integral and Pseudodifferential Equations with Numerical Approximation, Springer Monographs in Mathematics, SpringerVerlag, Berlin, Heidelberg, 2002.
 [4] G. W. Milton, The Theory of Composites, Vol. 6 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, UK, 2002.
 [5] V. Šmilauer, Z. Bittnar, Microstructurebased micromechanical prediction of elastic properties in hydrating cement paste, Cement and Concrete Research 36 (9) (2006) 1708–1718.
 [6] A. Prakash, R. A. Lebensohn, Simulation of micromechanical behavior of polycrystals: finite elements versus fast Fourier transforms, Modelling and Simulation in Materials Science and Engineering 17 (6) (2009) 064010+.
 [7] R. Brenner, J. BravoCastillero, Response of multiferroic composites inferred from a fastFouriertransformbased numerical scheme, Smart Materials and Structures 19 (11) (2010) 115004+.
 [8] Y. Saad, Iterative Methods for Sparse Linear Systems, second edition with corrections Edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2003.
 [9] J. Vondřejc, Analysis of heterogeneous materials using efficient meshless algorithms: Onedimensional study, Diploma thesis, Czech Technical University in Prague, URL: http://mech.fsv.cvut.cz/~vondrejc/download/ING.pdf (2009).