Green–Jacobi Preconditioner in FFT-Based Elasticity

• Green–Jacobi Preconditioner is a composite, symmetric approach that combines local Jacobi scaling with a global discrete Green operator to solve FFT-based periodic cell problems.
• It enhances convergence by mitigating the spectral dispersion issues of standard Green preconditioning in high-contrast, smoothly varying elasticity problems.
• The method preserves quasilinear computational complexity (O(N log N)) and adapts efficiently to applications like phase-field topology optimization and composite material simulations.

The Green–Jacobi preconditioner is a symmetric composite preconditioner for FFT-based solution of periodic cell problems on regular grids, introduced for small-strain linear elasticity in heterogeneous media whose coefficients are both smooth and highly contrasted. It is defined by

$M^{-1}=J^{1/2}\Gamma J^{1/2},$

where $\Gamma=(K^{\mathrm{ref}})^{-1}$ is the discrete Green operator of a constant reference medium and $J^{1/2}=(\operatorname{diag}K)^{-1/2}$ is the Jacobi factor extracted from the diagonal of the heterogeneous stiffness matrix. The construction is motivated by a specific failure mode of the standard Green preconditioner: although it yields nearly mesh-independent convergence for sharp, piecewise-constant microstructures, its performance degrades markedly for smoothly varying coefficient fields with high total contrast. The resulting Jacobi-accelerated FFT solver, termed J-FFT, preserves the quasilinear complexity $\mathcal{O}(N\log N)$ of conventional FFT-based solvers while substantially reducing conjugate-gradient iteration counts in the target regime (Ladecký et al., 4 Aug 2025).

1. Periodic cell problems and FFT-based linear elasticity

The setting is a periodic cell

$Y=\prod_{\alpha=1}^d [0,l_\alpha],$

equipped with periodic boundary conditions. In the small-strain formulation, the strain field is decomposed as

$$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$

where $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$ is the average strain and $u\in V$ is a periodic fluctuation field with zero mean. Using Mandel notation, the constitutive relation is

$\sigma(x)=C(x):\varepsilon(x),$

and equilibrium is expressed by

$-\nabla^\top \sigma(x)=0.$

The weak form is

$\Gamma=(K^{\mathrm{ref}})^{-1}$0

After finite-element discretization on a regular grid with piecewise linear basis functions, the equilibrium equations become

$\Gamma=(K^{\mathrm{ref}})^{-1}$1

where $\Gamma=(K^{\mathrm{ref}})^{-1}$2 collects nodal displacement unknowns, $\Gamma=(K^{\mathrm{ref}})^{-1}$3 is the discretized macroscopic strain, $\Gamma=(K^{\mathrm{ref}})^{-1}$4 is the discrete symmetrized-gradient operator, and $\Gamma=(K^{\mathrm{ref}})^{-1}$5 is the diagonal quadrature-weight matrix. For linear elasticity, Newton linearization yields

$\Gamma=(K^{\mathrm{ref}})^{-1}$6

In the numerical experiments, the stiffness tensor is taken as

$\Gamma=(K^{\mathrm{ref}})^{-1}$7

with fixed isotropic elasticity tensor $\Gamma=(K^{\mathrm{ref}})^{-1}$8 and scalar density parameter $\Gamma=(K^{\mathrm{ref}})^{-1}$9 controlling stiffness contrast (Ladecký et al., 4 Aug 2025).

This regular-grid discretization is central because it enables FFT acceleration. Under periodic boundary conditions and a constant reference medium, the corresponding discrete operator becomes block-circulant, and its Fourier representation is block-diagonal. That structure is the basis for both the standard Green preconditioner and its Green–Jacobi extension.

2. Standard discrete Green preconditioning and its failure mode

The standard FFT-based approach introduces a uniform reference medium through

$J^{1/2}=(\operatorname{diag}K)^{-1/2}$0

with spatially constant $J^{1/2}=(\operatorname{diag}K)^{-1/2}$1. Its inverse or pseudo-inverse,

$J^{1/2}=(\operatorname{diag}K)^{-1/2}$2

is the discrete Green operator. Because $J^{1/2}=(\operatorname{diag}K)^{-1/2}$3 is block-circulant on the periodic regular grid, $J^{1/2}=(\operatorname{diag}K)^{-1/2}$4 is applied as

$J^{1/2}=(\operatorname{diag}K)^{-1/2}$5

with $J^{1/2}=(\operatorname{diag}K)^{-1/2}$6 and $J^{1/2}=(\operatorname{diag}K)^{-1/2}$7 denoting the forward and inverse FFT and $J^{1/2}=(\operatorname{diag}K)^{-1/2}$8 block-diagonal in Fourier space. The standard preconditioned system is

$J^{1/2}=(\operatorname{diag}K)^{-1/2}$9

A single application of $\mathcal{O}(N\log N)$0 requires one FFT, one inverse FFT, and pointwise multiplication by $\mathcal{O}(N\log N)$1, so the complexity remains $\mathcal{O}(N\log N)$2 (Ladecký et al., 4 Aug 2025).

The paper identifies a precise limitation of this conventional construction. For discontinuous, piecewise-constant microstructures with relatively few phases, Green-preconditioned PCG exhibits nearly mesh-independent convergence, which is one reason FFT-based solvers are attractive on fine regular grids. However, when coefficients vary smoothly over the cell and the total contrast remains high, the iteration count can increase dramatically. The motivating example begins from a two-phase inclusion/matrix microstructure with contrast $\mathcal{O}(N\log N)$3, repeatedly smoothed by Gaussian filtering. For the initial sharp geometry, Green-PCG is efficient; as the field becomes smoother while overall contrast remains close to $\mathcal{O}(N\log N)$4, the Green-PCG iteration count increases strongly; only once smoothing becomes extreme enough that the total contrast drops substantially does the iteration count decrease again (Ladecký et al., 4 Aug 2025).

The paper emphasizes that the problematic regime is therefore not high contrast alone, but the combination of high total contrast and smoothly varying coefficients / many intermediate phases. Across laminate and cosine examples, two trends recur: iteration counts are almost independent of discretization size $\mathcal{O}(N\log N)$5, yet they increase with the number of material phases $\mathcal{O}(N\log N)$6, or more broadly with smoothness and multiplicity of coefficient values. A counterintuitive point follows: smooth data can be worse than sharp interfaces from the iterative-solver point of view. The discussion explicitly distinguishes this from discretization error; the issue concerns iterative/algebraic error, not the accuracy of the underlying spatial approximation (Ladecký et al., 4 Aug 2025).

The spectral rationale given is qualitative rather than a closed-form condition-number formula. Prior work cited in the paper indicates that Green/Laplace preconditioning yields a matrix close to diagonal, with eigenvalues related to local material properties. If the coefficient field spans many smoothly varying local values over a wide range, the preconditioned operator retains substantial eigenvalue spread. By contrast, a sharp two-phase configuration effectively associates the spectrum with only a small number of local coefficient values, which is more favorable for PCG.

3. Definition and interpretation of the Green–Jacobi preconditioner

The Jacobi preconditioner alone is

$\mathcal{O}(N\log N)$7

or explicitly

$\mathcal{O}(N\log N)$8

Used by itself, it provides only local scaling. It is cheap to apply, with $\mathcal{O}(N\log N)$9 cost, and it captures the heterogeneous diagonal amplitudes of the true stiffness matrix, but it lacks the global coupling information encoded by $Y=\prod_{\alpha=1}^d [0,l_\alpha],$0. Numerically, pure Jacobi is much worse than Green on fine grids because its iteration counts grow strongly with mesh size (Ladecký et al., 4 Aug 2025).

The Green–Jacobi construction combines these complementary effects. To preserve symmetry for PCG, the Jacobi factor is split as

$Y=\prod_{\alpha=1}^d [0,l_\alpha],$1

with

$Y=\prod_{\alpha=1}^d [0,l_\alpha],$2

The composite preconditioner is then defined by

$Y=\prod_{\alpha=1}^d [0,l_\alpha],$3

leading to the preconditioned system

$Y=\prod_{\alpha=1}^d [0,l_\alpha],$4

Using the Fourier representation of $Y=\prod_{\alpha=1}^d [0,l_\alpha],$5, the FFT-accelerated form becomes

$Y=\prod_{\alpha=1}^d [0,l_\alpha],$6

The paper terms the resulting method the Jacobi-accelerated FFT solver (J-FFT) (Ladecký et al., 4 Aug 2025).

The interpretation is straightforward. Green handles the global elliptic structure and yields mesh-scalable behavior; Jacobi corrects local coefficient scaling. In operational terms, Green–Jacobi applies a local rescaling in real space via $Y=\prod_{\alpha=1}^d [0,l_\alpha],$7, a global elliptic correction via $Y=\prod_{\alpha=1}^d [0,l_\alpha],$8 in Fourier space, and a second local rescaling via $Y=\prod_{\alpha=1}^d [0,l_\alpha],$9. From a linear-algebra viewpoint, PCG sees the symmetrically preconditioned operator

$$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$0

with $$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$1. The intended effect is tighter spectral clustering than for $$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$2, especially when stiffness values vary smoothly across a large range.

4. Algorithmic structure, complexity, and matrix-free assembly

The paper does not present full pseudocode, but it specifies the operator application. For a vector $$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$3,

$$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$4

is computed by five steps:

$$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$5

Then $$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$6. Matrix-vector products with the stiffness operator remain matrix-free:

$$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$7

The per-application cost is dominated by the FFT pair. Two diagonal multiplications by $$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$8 cost $$\varepsilon(x)=\overline{\varepsilon} + \nabla^s u(x),$$9, one FFT and one inverse FFT cost $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$0, and the block-diagonal multiplication by $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$1 costs $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$2. The resulting complexity is therefore

$$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$3

which preserves the quasilinear scaling that motivates FFT-based solvers in the first place (Ladecký et al., 4 Aug 2025).

A practical issue is that FFT-based FE solvers often use the operator $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$4 without assembling the full matrix $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$5. In that setting, $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$6 is not immediately available. A single diagonal entry can be extracted from

$$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$7

where $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$8 is a unit vector associated with degree of freedom $$\overline{\varepsilon}\in\mathbb{R}^{d_\varepsilon}$$9. A naive approach would require $u\in V$0 operator applications, which would be quadratic in cost. The paper instead exploits the sparsity and locality of the FE operator on the regular grid: for linear finite elements, the full diagonal can be assembled in

$u\in V$1

matrix-vector products, preserving overall linear complexity $u\in V$2 for Jacobi assembly. If some diagonal entries vanish, as may occur in the presence of voids, the recommended treatment is to replace zero diagonal entries by ones before inversion in order to avoid division by zero (Ladecký et al., 4 Aug 2025).

The implementation assumptions behind FFT acceleration remain unchanged: periodic boundary conditions, regular discretization, and a constant reference medium. The discrete Green operator is stored through its Fourier-space symbol $u\in V$3 rather than as a full matrix, while the additional storage cost of Green–Jacobi is only the diagonal of $u\in V$4, requiring $u\in V$5 memory.

5. Numerical behavior and comparative performance

The reported numerical evidence consistently separates three regimes: pure Green, pure Jacobi, and Green–Jacobi. Green remains nearly mesh-independent but deteriorates with smoothness and multiplicity of coefficient values. Jacobi is comparatively insensitive to the number of phases but scales poorly with mesh refinement. Green–Jacobi occupies an intermediate position: it reintroduces some mesh dependence, though much milder than pure Jacobi, while greatly reducing sensitivity to smooth heterogeneous coefficients (Ladecký et al., 4 Aug 2025).

In the Gaussian-filtered inclusion example at contrast $u\in V$6, Green-PCG requires fewer iterations than Green–Jacobi-PCG for the original sharp geometry. As filtering proceeds and the field becomes smoother while contrast remains high, Green-PCG iteration counts increase dramatically, whereas Green–Jacobi overtakes Green and remains efficient. The worst case for Green occurs at an intermediate level of smoothing, when the field is smooth but total contrast is still near $u\in V$7.

The laminate experiment uses

$u\in V$8

sampled with different numbers of material pixels $u\in V$9 and solved on FE grids $\sigma(x)=C(x):\varepsilon(x),$0. For Green preconditioning, iteration counts are stable with respect to $\sigma(x)=C(x):\varepsilon(x),$1 but grow with $\sigma(x)=C(x):\varepsilon(x),$2; at contrast $\sigma(x)=C(x):\varepsilon(x),$3 this growth saturates, whereas at contrast $\sigma(x)=C(x):\varepsilon(x),$4 it continues. For Jacobi, iteration counts grow strongly with $\sigma(x)=C(x):\varepsilon(x),$5 and only mildly with $\sigma(x)=C(x):\varepsilon(x),$6. For Green–Jacobi, there is some growth with $\sigma(x)=C(x):\varepsilon(x),$7, but much weaker than for Jacobi, and as $\sigma(x)=C(x):\varepsilon(x),$8 increases the iteration counts tend to decrease; even at contrast $\sigma(x)=C(x):\varepsilon(x),$9, the increase over contrast $-\nabla^\top \sigma(x)=0.$0 is less than a factor of two (Ladecký et al., 4 Aug 2025).

The cosine example,

$-\nabla^\top \sigma(x)=0.$1

is studied for $-\nabla^\top \sigma(x)=0.$2 and $-\nabla^\top \sigma(x)=0.$3, the latter introducing true voids. Green iteration counts again increase with the number of phases $-\nabla^\top \sigma(x)=0.$4, and slightly with $-\nabla^\top \sigma(x)=0.$5. Jacobi again shows poor mesh scaling. Green–Jacobi iteration counts increase with $-\nabla^\top \sigma(x)=0.$6, but much more slowly than Jacobi, and decrease monotonically as $-\nabla^\top \sigma(x)=0.$7 increases and the material field becomes smoother. The method also works in the presence of voids, indicating that its usefulness extends beyond merely large finite contrasts (Ladecký et al., 4 Aug 2025).

A practically important experiment concerns phase-field topology optimization, where the density field evolves from random noise to a near-binary microstructure with diffuse interfaces, and the equilibrium problem

$-\nabla^\top \sigma(x)=0.$8

must be solved repeatedly for several load cases. Over the L-BFGS optimization process on meshes $-\nabla^\top \sigma(x)=0.$9, $\Gamma=(K^{\mathrm{ref}})^{-1}$00, and $\Gamma=(K^{\mathrm{ref}})^{-1}$01, the Green preconditioner requires more than $\Gamma=(K^{\mathrm{ref}})^{-1}$02 iterations for $\Gamma=(K^{\mathrm{ref}})^{-1}$03, more than $\Gamma=(K^{\mathrm{ref}})^{-1}$04 for $\Gamma=(K^{\mathrm{ref}})^{-1}$05, and more than $\Gamma=(K^{\mathrm{ref}})^{-1}$06 for $\Gamma=(K^{\mathrm{ref}})^{-1}$07 late in optimization. Green–Jacobi requires approximately $\Gamma=(K^{\mathrm{ref}})^{-1}$08 iterations for $\Gamma=(K^{\mathrm{ref}})^{-1}$09, $\Gamma=(K^{\mathrm{ref}})^{-1}$10 for $\Gamma=(K^{\mathrm{ref}})^{-1}$11, and $\Gamma=(K^{\mathrm{ref}})^{-1}$12 for $\Gamma=(K^{\mathrm{ref}})^{-1}$13. Pure Jacobi also outperforms Green here, with less than $\Gamma=(K^{\mathrm{ref}})^{-1}$14 iterations for $\Gamma=(K^{\mathrm{ref}})^{-1}$15, less than $\Gamma=(K^{\mathrm{ref}})^{-1}$16 for $\Gamma=(K^{\mathrm{ref}})^{-1}$17, and less than $\Gamma=(K^{\mathrm{ref}})^{-1}$18 for $\Gamma=(K^{\mathrm{ref}})^{-1}$19, but Green–Jacobi remains clearly best among the three (Ladecký et al., 4 Aug 2025).

A final comparison isolates the role of interphase smoothness. For a topology-optimized design represented once as a smooth density field $\Gamma=(K^{\mathrm{ref}})^{-1}$20 and once as a thresholded sharp two-phase field $\Gamma=(K^{\mathrm{ref}})^{-1}$21, both solved at contrasts $\Gamma=(K^{\mathrm{ref}})^{-1}$22, $\Gamma=(K^{\mathrm{ref}})^{-1}$23, and $\Gamma=(K^{\mathrm{ref}})^{-1}$24, the smooth case exhibits slower convergence as contrast increases, and Green–Jacobi outperforms Green more strongly at higher contrast. In the sharp case, convergence is almost independent of total contrast, and Green outperforms Green–Jacobi. This identifies smoothness of interphases, rather than contrast alone, as the main reason the standard Green preconditioner loses efficiency.

6. Relevance, stopping criterion, and limitations

The paper emphasizes three application domains. In phase-field fracture, the scalar damage or phase field degrades stiffness smoothly over diffuse interfaces, creating large contrast together with smooth transitions. In density-based topology optimization, especially with phase-field regularization, material properties vary smoothly between void-like and solid-like states, and the elasticity problem must be solved repeatedly for changing coefficient fields. In adaptive-grid and composite-voxel FFT solvers, as well as XFEM-like enriched regular-grid schemes, sharp voxel jumps are often replaced by smoothed or mixed effective coefficients, increasing the number of local material states while maintaining substantial overall contrast. These are precisely the regimes in which Green–Jacobi is beneficial (Ladecký et al., 4 Aug 2025).

For solver comparison, the stopping criterion is based on the $\Gamma=(K^{\mathrm{ref}})^{-1}$25-norm of the residual,

$\Gamma=(K^{\mathrm{ref}})^{-1}$26

where the $\Gamma=(K^{\mathrm{ref}})^{-1}$27-inner product is

$\Gamma=(K^{\mathrm{ref}})^{-1}$28

The paper states that this criterion is preferred because it is related to energy error and error in homogenized properties.

The strengths and limitations are explicitly balanced. Green–Jacobi preserves FFT-based $\Gamma=(K^{\mathrm{ref}})^{-1}$29 complexity, significantly improves convergence for smooth high-contrast material data, is effective in practical applications such as phase-field topology optimization, works with voids or infinite contrast, and is straightforward to add to existing Green-preconditioned solvers. Its tradeoffs are equally clear: it is not as mesh-independent as pure Green preconditioning; for sharp two-phase microstructures with few phases, standard Green preconditioning can still be better; there is extra setup cost for assembling the Jacobi diagonal; and each iteration is slightly more expensive than Green-PCG because of the two additional diagonal scalings (Ladecký et al., 4 Aug 2025).

A plausible implication is that the Green–Jacobi preconditioner marks a shift in how FFT-based solvers are tuned for modern coefficient fields. Classical discrete Green preconditioning remains highly effective for binary or few-phase media with sharp interfaces, but the increasing use of diffuse interfaces, filtering, regularization, and mixed-cell representations creates a distinct algebraic regime. In that regime, local coefficient normalization through the Jacobi factor complements the global Fourier-space correction of the discrete Green operator, yielding a preconditioner specifically adapted to smooth, highly heterogeneous coefficients.