Analysis of a Fast Fourier Transform Based Method for Modeling of Heterogeneous Materials

Analysis of a Fast Fourier Transform Based Method for Modeling of Heterogeneous Materials

Jaroslav Vondřejc 1Czech Technical University in Prague,
Faculty of Civil Engineering, Department of Mechanics
   Jan Zeman 1Czech Technical University in Prague,
Faculty of Civil Engineering, Department of Mechanics
   Ivo Marek 2Czech Technical University in Prague,
Faculty of Civil Engineering, Department of Mathematics2

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 [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 Gradient-based scheme over the original Moulinec-Suquet algorithm.

Homogenization, Fast Fourier Transform, Conjugate Gradients

1 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 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 [1]. 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 [2], we have offered a new approach to the Moulinec-Suquet scheme, by exploiting the trigonometric collocation method due to Saranen and Vainikko [3]. 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.

The original problem (2)–(3) is then equivalent to the periodic Lippmann-Schwinger integral equation, formally written as


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


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 [5], (visco-)plasticity [6] or to multiferroics [7].

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. [2]


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 [1] 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 [2] 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 .

Definition 1

Given , we define operator and associated sub-space as

Lemma 1

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 [8]. ∎

Remark 1

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].

Lemma 2

The solution to the linear system (9) admits the decomposition , with satisfying


As , Lemma 1 ensures that it can be decomposed into two orthogonal parts and . Substituting this expression into (9), and using the identity , we arrive at


Since , we have and the proof is complete. ∎

With these auxiliary results in hand, we are in the position to present our main result.

Proposition 1

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 .

Proof (outline)

It follows from Lemma 2 that the solution to (9) admits yet another, optimization-based, characterization in the form


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 [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 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 [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 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.

Figure 1: (a) Relative number of iterations as a function of the reference medium parameter and (b) computational time as a function of the number of unknowns.

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  [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 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.

5 Conclusions

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:

  1. 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.

  2. 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.


  • [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 FFT-based 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, Springer-Verlag, 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, Microstructure-based 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. Bravo-Castillero, Response of multiferroic composites inferred from a fast-Fourier-transform-based 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: One-dimensional study, Diploma thesis, Czech Technical University in Prague, URL: (2009).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description