An FFT-based Galerkin Method for Homogenization of Periodic Media

Jaroslav Vondřejc; Jan Zeman; Ivo Marek

Source-linked AI summary

An FFT-based Galerkin Method for Homogenization of Periodic Media

Jaroslav Vondřejc, Jan Zeman, Ivo Marek

arXiv:1311.0089v2Nov 1, 2013math.NAcond-mat.mtrl-scimath.APphysics.comp-ph

TL;DR

The paper addresses how to understand and improve FFT-based solvers for periodic homogenization cell problems. It recasts the Moulinec–Suquet scheme as a trigonometric-polynomial Galerkin discretization, proves convergence in the scalar elliptic setting, and explains why conjugate gradients applies despite a non-symmetric system.

Problem
Finite element solution of periodic cell problems can become computationally prohibitive for high-resolution coefficient data, motivating analysis of efficient FFT-based methods.
Method
The paper formulates the cell problem variationally, discretizes it with trigonometric-polynomial Galerkin spaces and numerical integration, and analyzes the resulting linear system.
Results
The discretized solutions converge to the weak solution, and the numerically integrated Galerkin method is equivalent to the discrete Lippmann–Schwinger equation underlying Moulinec–Suquet.
Takeaways & Limitations
The Galerkin interpretation provides theoretical support for FFT-based numerical homogenization and explains the applicability of conjugate gradients to the resulting system.
Takeaways & Limitations
The analysis assumes uniformly elliptic, symmetric coefficients and focuses on the scalar elliptic setting.

Abstract

from arXiv · show

In 1994, Moulinec and Suquet introduced an efficient technique for the numerical resolution of the cell problem arising in homogenization of periodic media. The scheme is based on a fixed-point iterative solution to an integral equation of the Lippmann-Schwinger type, with action of its kernel efficiently evaluated by the Fast Fourier Transform techniques. The aim of this work is to demonstrate that the Moulinec-Suquet setting is actually equivalent to a Galerkin discretization of the cell problem, based on approximation spaces spanned by trigonometric polynomials and a suitable numerical integration scheme. For the latter framework and scalar elliptic setting, we prove convergence of the approximate solution to the weak solution, including a-priori estimates for the rate of convergence for sufficiently regular data and the effects of numerical integration. Moreover, we also show that the variational structure implies that the resulting non-symmetric system of linear equations can be solved by the conjugate gradient method. Apart from providing a theoretical support to Fast Fourier Transform-based methods for numerical homogenization, these findings significantly improve on the performance of the original solver and pave the way to similar developments for its many generalizations proposed in the literature.

1. Introduction

The paper frames FFT-based homogenization as an efficient alternative for periodic cell problems that can become costly for high-resolution material data. It explains the method’s convergence behavior and establishes its equivalence to a trigonometric-polynomial Galerkin discretization.

Problem setting: Periodic homogenization requires solving a cell problem involving periodic gradient, flux, and material-coefficient fields.The gradient field combines a mean applied gradient with a fluctuating periodic component.
Motivation: Finite element treatment can become computationally prohibitive when coefficients come from large, high-resolution imaging datasets.The paper motivates FFT-based methods specifically for such coefficient data.
Existing FFT solver: Moulinec–Suquet reformulates the cell problem as a Lippmann–Schwinger equation and evaluates its Green-operator action with FFTs.The numerical scheme uses fixed-point iterations.
Existing FFT solver: The original solver’s convergence depends on the reference tensor, while iteration counts grow linearly with coefficient contrast.Later work identified the reference-tensor dependence and contrast scaling.
Paper contribution: The authors explain earlier observations by identifying the Moulinec–Suquet discretization with a Galerkin method using trigonometric-polynomial approximation spaces.The paper also studies convergence, numerical integration, and solver properties in the scalar elliptic setting.
Paper contribution: The resulting non-symmetric linear system can nevertheless be solved by the standard conjugate-gradient algorithm.This observation motivates the paper’s variational explanation of the solver behavior.

2. Notation and preliminaries

This section establishes the notation, periodic function spaces, tensor conventions, and Fourier tools used throughout the Galerkin analysis. Trigonometric and Sobolev-space machinery provides the framework for representing periodic fields and their regularity.

Role in the paper: The section’s notation and Fourier-transform facts prepare the rigorous formulation of the Lippmann–Schwinger equation and weak cell problem.These tools support the projection-based equivalence developed later.
Algebraic notation: Vectors, tensors, inner products, tensor products, repeated-index summation, and the Kronecker delta are defined for compact formulas.The unit tensor is represented using the Kronecker delta.
Function-space notation: The paper introduces notation for periodic scalar-, vector-, and tensor-valued functions and their standard norms.The notation includes Lp spaces, essential suprema, and the measure of the periodic cell.
Fourier representation: Periodic Fourier basis functions are introduced as an orthonormal basis for the relevant L2 spaces.Their wave-vector notation is used in subsequent Fourier-domain arguments.
Regularity spaces: Periodic Sobolev spaces Hs are defined with Fourier-based norms, alongside spaces of differentiable periodic functions.Sobolev inequalities are invoked for sufficiently large regularity indices.

3. Weak and integral formulations of cell problem

The paper rigorously connects the weak cell problem with its Lippmann–Schwinger formulation through a projection operator. Under ellipticity and symmetry assumptions, the formulations are equivalent and uniquely solvable, enabling Galerkin discretization.

3.1. Problem setting: The weak and integral formulations are developed in periodic spaces using divergence, curl, bilinear forms, and a projection operator.The projection reflects the differential constraints of the cell problem.
3.1. Problem setting: Periodic vector fields admit an orthogonal decomposition into constant, curl-free, and divergence-free components.This Helmholtz-type decomposition organizes the constrained solution space.
3.1. Problem setting: The coefficient tensor is assumed essentially bounded, symmetric, and uniformly elliptic, with contrast measured by ρ_A = C_A/c_A.The auxiliary tensor A^(0) is subject to corresponding structural assumptions.
3.2. Projection operator: The projection properties follow from Fourier representations, symmetry of A^(0), and the Helmholtz decomposition.These properties establish the operator’s boundedness, invariance, and orthogonality in the scalar reference case.
3.3. Equivalence of solutions: The weak cell problem and Lippmann–Schwinger equation have equivalent unique solutions under assumptions (A1)–(A3).The equivalence identifies the fluctuating field with the projected component and the mean field with the complementary component.

4. Discretization

The paper discretizes the cell problem by Galerkin projection onto trigonometric-polynomial spaces, whose Fourier and grid representations support FFT computation. It establishes qualitative and convergence results for exact and numerically integrated formulations under stated regularity and grid assumptions.

4.1. Trigonometric polynomials: The discrete representation uses a regular N1 × ... × Nd grid, grid spacings hα, and the Discrete Fourier Transform to connect grid and Fourier values.The grid is assumed symmetric about the origin, with every Nα odd.
4.1. Trigonometric polynomials: The truncation operator PN is orthogonal, whereas the interpolation operator QN is a projection but not an orthogonal one.Their approximation properties underpin the trigonometric-polynomial discretization.
4.1. Trigonometric polynomials: Trigonometric polynomials provide conforming finite-dimensional approximations of the relevant constant, curl-free, and divergence-free spaces.The construction preserves the structure of the continuous spaces and supports a trigonometric Helmholtz decomposition.
4.2. Galerkin approximation: Galerkin projection of the variational cell problem has a unique solution under the stated coefficient and grid assumptions, with error estimates obtained through the Cea lemma.The convergence argument uses approximation estimates for trigonometric-polynomial projections.
4.3. Galerkin approximation with numerical integration: The numerical-integration discretization relies on coefficient and data conditions ensuring that the interpolated bilinear and linear forms are well-defined.The paper notes that exact form evaluation is generally unavailable for arbitrary coefficients, motivating numerical integration.
4.3. Galerkin approximation with numerical integration: With numerical integration, the discrete problem remains uniquely solvable under additional assumptions, while convergence-rate estimates require sufficient regularity of the weak solution.The numerical-integration analysis uses the interpolation operator QN and Strang-lemma estimates.

5. Algebraic system and its solution

The fully discrete GaNi formulation is represented with structured grid and Fourier-domain objects, and is equivalent to the discrete Lippmann–Schwinger system. Its variational structure makes the resulting non-symmetric system solvable by conjugate gradients, with FFT-based matrix actions.

5.1. Notation and preliminaries: Structured vectors and matrices encode grid values and their Fourier representations for the fully discrete formulation.The notation includes spaces for constant, zero-mean curl-free, and divergence-free trigonometric-polynomial fields.
5.2. Fully discrete formulations: The interpolation-based numerical integration framework converts the variational problem into a fully discrete GaNi system.The operator IN connects trigonometric polynomials with their grid values and is one-to-one and isometric under the stated assumptions.
5.2. Fully discrete formulations: The discrete GaNi solution, discrete Lippmann–Schwinger solution, and equivalent linear system coincide under the stated assumptions.For A(0) = λI with λ > 0, the discrete variational formulation yields the corresponding linear system.
5.3. Solution of linear system: FFT evaluation gives the structured system an action cost of O(|N| log |N|), enabling efficient iterative solution.The matrix is non-symmetric but consists of sparse structured factors whose dominant operations are forward and inverse discrete Fourier transforms.
5.3. Solution of linear system: The non-symmetric system can be solved by conjugate gradients for arbitrary initial solutions in the discrete curl-free space.Self-adjointness and positive-definiteness on the relevant subspace keep all iterates within that space.
5.3. Solution of linear system: Conjugate-gradient iterations for a fixed tolerance grow as √ρA, whereas the original fixed-point scheme can grow linearly with coefficient contrast.The conjugate-gradient result follows from κ(AN) = ρA; the comparison concerns the original scheme's contrast dependence.

6. Conclusions

The paper establishes a trigonometric-polynomial Galerkin framework for the scalar elliptic cell problem and connects it to the Moulinec–Suquet discretization. The analysis proves convergence, identifies the discrete equivalence, and enables conjugate-gradient solution of the resulting system.

6. Conclusions: Trigonometric polynomials provide conforming, structure-preserving approximations to infinite-dimensional curl-free spaces.The approximation uses a finite-dimensional space and a projection operator reflecting the problem's differential constraints.
6. Conclusions: Solutions of the discretized problems with or without numerical integration converge to the weak solution at standard rates for sufficiently regular data.The conclusion covers both discretization variants and regularity-dependent convergence rates.
6. Conclusions: The GaNi method is equivalent to the discrete Lippmann–Schwinger equation underlying the original Moulinec–Suquet scheme.This equivalence supplies a variational interpretation of the FFT-based discretization.
6. Conclusions: The non-symmetric GaNi linear system is independent of the auxiliary parameter A(0) and can be solved by conjugate gradients.These properties are stated as consequences of the variational formulation and its discrete structure.
6. Conclusions: The framework is presented as a starting point for extensions to error estimation, solver refinements, and more complex physical or material settings.The paper describes these directions as possibilities for future investigations rather than established results.

7. Comparison with results by Brisard and Dormieux [9]

The paper contrasts its curl-free trigonometric-polynomial Galerkin framework with the Brisard–Dormieux approach based on polarization fields and pixel- or voxel-wise constants. The schemes differ in consistency, conformity, parameter dependence, and iterative-solver requirements.

7. Comparison with results by Brisard and Dormieux [9]: Brisard and Dormieux formulate the discretization through stationarity conditions for a Hashin–Shtrikman functional in an unknown polarization field.Their polarization field lies in the full L2 space and is related to the Lippmann–Schwinger formulation.
7. Comparison with results by Brisard and Dormieux [9]: Their approximation uses pixel- or voxel-wise constant fields, allowing the formulation to be localized to individual pixels or voxels.This differs from the present work's approximation in the zero-mean curl-free subspace.
7. Comparison with results by Brisard and Dormieux [9]: Both approaches consider consistent and non-consistent approximations, but the earlier work does not provide a-priori convergence-rate estimates.The non-consistent approximation truncates the relevant infinite sum, with an effect comparable to numerical integration.
7. Comparison with results by Brisard and Dormieux [9]: The consistent Green operator requires highly accurate lattice-sum evaluation, which is difficult especially in three dimensions.Truncating the lattice sum instead introduces numerical-integration errors.
7. Comparison with results by Brisard and Dormieux [9]: Truncation produces non-conforming gradient fields, so the discrete Helmholtz decomposition no longer holds.The loss of conformity is identified as a consequence of truncating the lattice sum.
7. Comparison with results by Brisard and Dormieux [9]: In the earlier formulation, the system matrix depends on A(0) and may be positive-definite, negative-definite, or indefinite.The stationary point's classification depends on the auxiliary parameter, requiring more complex iterative solvers.

Appendix A. Approximation by trigonometric polynomials

Appendix A establishes approximation and operator representations for trigonometric-polynomial spaces, including convergence conditions for a key Fourier sum. It derives these properties using density, projection structure, and discrete Fourier relations.

Generalization: The results generalize established one-dimensional estimates to multidimensional vector settings with distinct grid spacings.The section explicitly frames these as generalizations of results from.
Approximation: Convergence in (25) follows from the density of trigonometric polynomials in L2.The proof invokes the density of {ϕk}k∈Zd in L2.
Operator representation: The operator QN is represented in Fourier space by using its projection property on trigonometric-polynomial spaces and element-by-element multiplication.The derivation proceeds through Fourier coefficients, linearity, and the projection structure of QN.
Estimates: The estimates combine previously established bounds, difference-quotient relations, and the Cauchy inequality.The proof also uses a multidimensional extension of a one-dimensional estimate.

Appendix B. Regularity result

Appendix B develops periodic regularity results needed for convergence-rate estimates of Galerkin approximations. Using difference quotients and periodic Sobolev arguments, it establishes H1 regularity and associated a-priori bounds for the weak cell-problem solution.

Purpose and setting: The appendix collects regularity results needed to establish Galerkin convergence rates for the periodic cell problem.The treatment uses difference-quotient techniques simplified by periodicity and assumes A ∈ W 1,∞.
Difference quotients: Difference quotients are controlled uniformly in the increment and coordinate index, enabling periodic Sobolev regularity arguments.The appendix states that suitable constants are independent of △ and α.
Regularity: The weak solution to the cell problem belongs to H1.This regularity follows after bounding the difference quotient independently of △.
Proof strategy: The regularity argument relies on translated coefficients, integration by parts for difference quotients, and induction for higher-order conclusions.The coefficient translation relation and induction-based extension are stated in the appendix.
A-priori estimates: Standard a-priori estimates yield inequality (B.1) once H1 regularity has been established.The stated estimate is obtained from standard bounds on the solution.