Analysis of a Fast Fourier Transform Based Method for Modeling of Heterogeneous Materials
The focus of this paper is on the analysis of the Conjugate Gradient method applied to a non-symmetric system of linear equations, arising from a Fast Fourier Transform-based homogenization method due to Moulinec and Suquet . 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 Gradient-based scheme over the original Moulinec-Suquet algorithm.
Keywords:Homogenization, Fast Fourier Transform, Conjugate Gradients
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 voxel-based geometry. Such data now allow for the design of bottom-up 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 large-scale voxel-based data in an efficient manner. In the engineering community, perhaps the most successful solver meeting these criteria was proposed by Moulinec and Suquet in . The algorithm is based on the Neumann series expansion of the inverse of an operator arising in the associated Lippmann-Schwinger equation and exploits the Fast Fourier Transform to evaluate the action of the operator efficiently for voxel-based data. In our recent work , we have offered a new approach to the Moulinec-Suquet scheme, by exploiting the trigonometric collocation method due to Saranen and Vainikko . Here, the Lippman-Schwinger equation is projected to a space of trigonometric polynomials to yield a non-symmetric 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 second-order 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 multi-index 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
Sub-matrices of are denoted as
for and . Analogously, the block vectors are denoted by lower case letters, e.g. and the matrix-by-vector multiplication is defined as
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
where is a -periodic vectorial electric field, denotes the corresponding vector of electric current and is a second-order positive-definite tensor of electric conductivity. In addition, the field is subject to a constraint
where denotes the -dimensional measure of and a prescribed macroscopic electric field.
Operator stands for the Fourier coefficient of for the -th frequency given by
”” 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 , (visco-)plasticity  or to multiferroics .
2.1 Discretization via trigonometric collocation
The numerical solution of the Lippmann-Schwinger 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])
where designates the Fourier coefficients defined in (6). Notice that the trigonometrical polynomials are uniquely determined by a regular grid data, which makes them well-suited to problems with pixel- or voxel-based computations.
The trigonometric collocation method is based on the projection of the Lippmann-Schwinger equation (4) onto the space of the trigonometric polynomials
leading to a to linear system in the form, cf. 
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 block-diagonal structure. In particular,
with , and . The matrix implements the Discrete Fourier Transform and is defined as
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) well-suited for applying some iterative solution technique. In particular, the original Fast Fourier Transform-based Homogenization scheme formulated by Moulinec and Suquet in  is based on the Neumann expansion of the matrix inverse , so as to yield the -th iterate in the form
As indicated earlier, our numerical experiments  suggest that the system can be efficiently solved using the Conjugate Gradient method, despite the non-symmetry 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 sub-space introduced later. Note that for simplicity, the reference conductivity is taken as .
Given , we define operator and associated sub-space as
The operator is an orthogonal projection.
First, we will prove that is projection, i.e. . Since is a unitary matrix, it is easy to see that
Hence, in view of the block-diagonal character of , it it sufficient to prove the projection property of sub-matrices only. This follows using a simple algebra, recall Eq. (5):
The orthogonality of now follows from
according to a well-known result of linear algebra, e.g. Proposition 1.8 in . ∎
It follows from the previous results that the subspace collects the non-zero 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].
The solution to the linear system (9) admits the decomposition , with satisfying
With these auxiliary results in hand, we are in the position to present our main result.
The non-symmetric 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 .
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 positive-definite on , the convergence of algorithm now follows from standard arguments, e.g. Theorem 6.6 in . Observe that different choices of generate identical Krylov subspaces, thus the sequence of iterates is independent of .∎
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 .
4 Numerical example
To support our theoretical results, we consider a three-dimensional model problem of electric conduction in a cubic periodic unit cell , representing a two-phase 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 nodes111In particular, was taken consequently as and leading up to unknowns. The conductivity of the homogeneous reference medium is parametrized as
where delivers the optimal convergence of the original Moulinec-Suquet Fast-Fourier Transform-based Homogenization () algorithm .
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 mildly-contrasted composite down to for . The minor differences visible especially for can be therefore attributed to accumulation of round-off errors. These observations fully confirm our theoretical results presented earlier in Section 3.
In Fig. 1, we present the total computational time222The problem was solved with a Matlab® in-house 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 . 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 short-recurrence CG-based scheme over alternative schemes for non-symmetric systems, such as GMRES. Finally note that finer discretizations can be treated by a straightforward parallel implementation.
In this work, we have proven the convergence of Conjugate Gradient method for a non-symmetric 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 Gradient-based 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.
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.
-  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.
-  J. Zeman, J. Vondřejc, J. Novák, I. Marek, Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients, Journal of Computational Physics 229 (21) (2010) 8065–8071. arXiv:1004.1122.
-  J. Saranen, G. Vainikko, Periodic Integral and Pseudodifferential Equations with Numerical Approximation, Springer Monographs in Mathematics, Springer-Verlag, Berlin, Heidelberg, 2002.
-  G. W. Milton, The Theory of Composites, Vol. 6 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, UK, 2002.
-  V. Šmilauer, Z. Bittnar, Microstructure-based micromechanical prediction of elastic properties in hydrating cement paste, Cement and Concrete Research 36 (9) (2006) 1708–1718.
-  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+.
-  R. Brenner, J. Bravo-Castillero, Response of multiferroic composites inferred from a fast-Fourier-transform-based numerical scheme, Smart Materials and Structures 19 (11) (2010) 115004+.
-  Y. Saad, Iterative Methods for Sparse Linear Systems, second edition with corrections Edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2003.
-  J. Vondřejc, Analysis of heterogeneous materials using efficient meshless algorithms: One-dimensional study, Diploma thesis, Czech Technical University in Prague, URL: http://mech.fsv.cvut.cz/~vondrejc/download/ING.pdf (2009).