<figure><img src="image_1.png" /></figure>

HAL
open science

# Accelerating the Adaptive Eyre-Milton FFT-based method for infinitely double contrasted media

Martin Dolbeau, Jérémy Bleyer, Karam Sab

► To cite this version:

Martin Dolbeau, Jérémy Bleyer, Karam Sab. Accelerating the Adaptive Eyre-Milton FFT-based method for infinitely double contrasted media. Comptes Rendus. Mécanique, 2024, 352 (G1), pp.251-267. 10.5802/crmeca.269. hal-04802224

HAL Id: hal-04802224

https://enpc.hal.science/hal-04802224v1

Submitted on 25 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Comptes Rendus
***
Mécanique

Martin Dolbeau, Jérémy Bleyer and Karam Sab

Accelerating the Adaptive Eyre-Milton FFT-based method for infinitely double
contrasted media

Volume 352 (2024), p. 251-267

Online since: 22 November 2024

https://doi.org/10.5802/crmeca.269

<figure><img src="image_3.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article / *Article de recherche*

# Accelerating the Adaptive Eyre-Milton FFT-based method for infinitely double contrasted media

## *Accélération du schéma Eyre-Milton adaptatif pour l'homogénéisation par FFT des milieux à double contraste infini*

Martin Dolbeau<sup>a</sup>, Jérémy Bleyer<sup>b</sup> and Karam Sab<sup>*,b</sup>

**Abstract.** Sab et al. (2024) have recently proposed an FFT-based iterative algorithm, termed Adaptive Eyre-Milton (AEM), for solving the Lippmann-Schwinger equation in the context of periodic homogenization of infinitely double contrasted linear elastic composites (heterogeneous materials with linear constitutive laws that contain both pores and rigid inclusions). They have demonstrated the unconditional linear convergence of this scheme, regardless of initialization and the chosen reference material. However, numerical simulations have shown that the rate of convergence of AEM strongly depends on the chosen reference material. In this paper, we introduce a new version of the AEM scheme where the reference material is updated iteratively, resulting in a fast and versatile scheme, termed Accelerated Adaptive Eyre-Milton (A2EM). Numerical simulations with A2EM on linear elastic composites with both pores and infinitely rigid inclusions show that, regardless of the initial chosen reference material, this algorithm has the same rate of convergence as AEM with the best choice of reference material.

**Résumé.** Sab et al. (2024) ont récemment proposé un algorithme itératif basé sur la FFT, appelé Adaptive Eyre-Milton (AEM), pour résoudre l'équation de Lippmann-Schwinger dans le contexte de l'homogénéisation périodique de composites élastiques linéaires à double contraste infini (matériaux hétérogènes avec des lois constitutives linéaires contenant à la fois des pores et des inclusions rigides). Ils ont démontré la convergence linéaire inconditionnelle de ce schéma, quel que soit l'initialisation et le matériau de référence choisi. Cependant, les simulations numériques ont montré que la vitesse de convergence du schéma AEM dépend fortement du choix du matériau de référence. Dans cet article, nous introduisons une nouvelle version du schéma AEM, où le matériau de référence est mis à jour de manière itérative, aboutissant à un schéma rapide et polyvalent, appelé Accelerated Adaptive Eyre-Milton (A2EM). Des simulations numériques avec le schéma A2EM sur des composites élastiques linéaires avec des pores et des inclusions infiniment rigides montrent que, quel que soit le matériau de référence initial choisi, cet algorithme a la même vitesse de convergence que le schéma AEM avec le meilleur choix du matériau de référence.

**Keywords.** computational homogenization, FFT-based method, iterative scheme, linear elasticity, composite materials.

**Mots-clés.** homogénéisation computationnelle, méthode basée sur la FFT, schéma itératif, élasticité linéaire, matériaux composites.

Manuscript received 11 July 2024, accepted 26 September 2024.

<footnote><sup>a</sup> Laboratoire Navier, IPParis ENPC, Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France</footnote>

<footnote>*E-mails*: martin.dolbeau@enpc.fr, jeremy.bleyer@enpc.fr, karam.sab@enpc.fr</footnote>

<footnote>*Corresponding author</footnote>
# 1. Introduction

Since Moulinec and Suquet [1, 2] pioneering work, FFT-based methods for periodic homogenization of composites have become very popular. The initial corrector problem on the unit cell is reformulated as a single integral equation known as the Lippmann-Schwinger equation [3-5], which is solved iteratively. At each iteration, the convolution kernel of the integral equation must be applied. This is done most efficiently in the Fourier space, by means of fast Fourier transforms (FFT) in a discrete setting. Many authors introduced new iterative schemes which prove more efficient in some situations. Among them, the Eyre and Milton scheme [6], the *augmented Lagrangian* scheme of Michel, Moulinec and Suquet [7], the *polarization-based* schemes of Monchiet and Bonnet [8, 9]. In a non-linear setting, Newton or quasi-Newton approaches, combined with Krylov solvers [10, 11] deliver very efficient solution schemes [12]. Most notably, by observing that the basic scheme could be seen as a gradient descent, Kabel et al. [13] allowed the introduction of accelerated-gradient methods and other well-known optimization methods [14, 15]. The literature on the topic is very rich, and the reader is referred to the recent and comprehensive review by Schneider [16].

However, existing FFT-based algorithms were not adapted to the case of infinitely double contrasted linear elastic composites (heterogeneous materials with linear constitutive laws that contain both pores and rigid inclusions). Recently, the authors introduced in [17] the so-called polarization-based *adaptive Eyre-Milton* (AEM) scheme to address this case. They demonstrate the unconditional linear convergence of this scheme, regardless of initialization and the chosen reference material. Their numerical experiments indeed confirmed that AEM convergence rate is very fast if the reference material is well chosen.

In this paper, numerical simulations with AEM, and another polarization-based scheme named “$\gamma = 0$”-scheme, on linear elastic composites consisting of a homogeneous matrix containing both pores and infinitely rigid inclusions will demonstrate the existence of an optimal reference material. Then, an *accelerated* version of the AEM and “$\gamma = 0$” schemes where the reference material is updated iteratively will be presented. Two possible updating rules will be considered and compared. Numerical simulations will show that the rate of convergence of both the termed Accelerated Adaptive Eyre-Milton (A2EM) scheme and accelerated “$\gamma = 0$”-scheme does not depend much on the chosen initial reference material, and that it is more or less equal to the rate of convergence of the non accelerated original scheme with the best choice of reference material.

Following [17], the homogenization problem in the presence of pores and rigid inclusions is formulated in Section 2, as well as the corresponding Lippmann-Schwinger equation. The AEM scheme is recalled in Section 3. Then, the knew A2EM scheme is introduced in Section 4. This section is followed by the discussion in Section 5 of a few practical issues regarding the implementation of the method, namely: the optimization of the number of convolutions per iteration and the number of variables to be stored.

The paper closes in Section 6 with a few illustrative applications, in a three-dimensional setting. These examples clearly illustrate the performance of the proposed A2EM scheme.

# 2. The homogenization problem

We consider homogenization of a periodic, linearly elastic material, containing pores and rigid inclusions, in the $d$-dimensional real space $\mathbb{R}^d(d = 2,3)$. The unit-cell is denoted by $\Omega = [0, L_1] \times \dots \times [0, L_d]$ where $L_i > 0$ are the side lengths of $\Omega$ and $\langle \cdot \rangle_{\Omega}$ is the volume average operator over $\Omega$. Let $\Omega^f, \Omega^p$ and $\Omega^m$ denote the open domains of $\Omega$ occupied by, respectively, the rigid inclusions, the pores and the heterogeneous matrix.
A field $X$ is $\Omega$-periodic if $X(y_1, ..., y_d) = X(y_1 + L_1, ..., y_d) = ... = X(y_1, ..., y_d + L_d)$ at any point $\mathbf{y} = (y_1, ..., y_d) \in \mathbb{R}^d$. Let $L_{\text{sym}}^2(\Omega)$ denote the space of symmetric second-order tensor fields $\mathbf{y} \mapsto \mathbf{t}(\mathbf{y})$ ($t_{ij}(\mathbf{y}) = t_{ji}(\mathbf{y})$, $i, j = 1, ..., d$), which are $\Omega$-periodic and square integrable over $\Omega$. Adopting Einstein's summation convention over repeated indices, we introduce the double contraction operator “:” defined as $\mathbf{a} : \mathbf{b} = a_{ij} b_{ij}$, where $\mathbf{a}$ and $\mathbf{b}$ are two second-order tensors. The bilinear form $(\mathbf{a}, \mathbf{b}) \mapsto \langle \mathbf{a} : \mathbf{b} \rangle_{\Omega}$ then defines a scalar product over $L_{\text{sym}}^2(\Omega)$ and the associated natural norm is

$$ \|\mathbf{t}\|_{L_{\text{sym}}^2(\Omega)} = \sqrt{\langle \mathbf{t} : \mathbf{t} \rangle_{\Omega}}, \qquad (1) $$

It is well-known (see e.g. [18]) that, endowed with this natural norm and scalar product, $L_{\text{sym}}^2(\Omega)$ can be decomposed into the two orthogonal sub-spaces $S$ and $D$:

$$ L_{\text{sym}}^2(\Omega) = S \oplus D \qquad (2) $$

where $S$ is the sub-space of the $\sigma \in L_{\text{sym}}^2(\Omega)$ which are divergence-free in the sense of distributions over $\mathbb{R}^d$, $\sigma_{ij,j} = 0$, and $D$ is the sub-space of the strains $\mathbf{e} \in L_{\text{sym}}^2(\Omega)$ which are kinematically compatible with zero mean ($\mathbf{e}$ is curl-curl-free and $e_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i})$ where $\mathbf{u} = (u_i)$ is an $\Omega$-periodic displacement vector field).

Within the framework of linear elasticity, the fourth-order stiffness tensor $C(\mathbf{y}) = (C_{ijkl}(\mathbf{y}))$ ($i, j, k, l = 1, ..., d$) is defined only in the heterogeneous matrix for $\mathbf{y} \in \Omega^m$, and it exhibits both minor ($C_{ijkl} = C_{ijkl} = C_{ijkl}$) and major ($C_{ijkl} = C_{ijkl}$) symmetries.

We say that a pair of strain-stress fields, $(\varepsilon, \sigma) \in L_{\text{sym}}^2(\Omega) \times L_{\text{sym}}^2(\Omega)$, complies with the constitutive equations if the following equations hold true:

$$ \varepsilon = 0 \text{ in } \Omega^r, \quad \sigma = 0 \text{ in } \Omega^p \quad \text{and} \quad \sigma = C : \varepsilon \text{ in } \Omega^m. \qquad (3) $$

Let $E$ denote a macroscopic strain tensor (symmetric, second-order tensor which is uniform over $\Omega$). Then, the homogenization problem can be stated as follows:

$$ \text{Find } (\varepsilon_E, \sigma_E) \in D \times S \text{ such that } (\varepsilon_E = E + \varepsilon_E, \sigma_E) \text{ complies with (3)} \qquad (4) $$

It is shown that under the assumptions stated in [17] the above problem has a unique solution, $(\varepsilon_E, \sigma_E)$, up to a pair $(\mathbf{e}, \sigma) \in D \times S$ with $(\mathbf{e}, \sigma) = (0, 0)$ in $\Omega^m$, and that the homogenized stiffness tensor $C_{\text{hom}}$ is defined by:

$$ \forall E, \frac{1}{2} E : C_{\text{hom}} : E = \left\langle \frac{1}{2} \sigma_E : e_E \right\rangle_{\Omega} \qquad (5) $$

It is also shown that the problem (4) can be equivalently reformulated as follows:

$$ \text{Find } (\varepsilon_E, \sigma_E) \in L_{\text{sym}}^2(\Omega) \times L_{\text{sym}}^2(\Omega) \text{ complying to (3) such that } X(\varepsilon_E, \sigma_E) = 0, \qquad (6) $$

where, for any pair of strain-stress fields $(\varepsilon, \sigma) \in L_{\text{sym}}^2(\Omega) \times L_{\text{sym}}^2(\Omega)$, the corresponding residual field $X(\varepsilon, \sigma)$ is given by:

$$ X(\varepsilon, \sigma) = E - \varepsilon - \Gamma^0 * [(\sigma - C_0 : \varepsilon)]. \qquad (7) $$

Here, $C_0$ is the reference material which is a uniform, positive-definite $(\varepsilon : C_0 : \varepsilon > 0$ for all $\varepsilon \neq 0)$, fourth-order tensor with both minor and major symmetries, and $\Gamma^0$ its associated Green operator defined as follows.

For any $\tau \in L_{\text{sym}}^2(\Omega)$, the following problem:

$$ \text{Find } \mathbf{e}_{\tau} \in D \text{ such that } C_0 : \mathbf{e}_{\tau} + \tau \in S \qquad (8) $$

has a unique solution $\mathbf{e}_{\tau}$ that depends linearly on $\tau$. The Green operator $\Gamma^0$ associated to $C_0$ is defined as the linear operator that maps $\tau$ onto $-\mathbf{e}_{\tau}$ [3–5]:

$$ \mathbf{e}_{\tau} = -\Gamma^0 * \tau \qquad (9) $$
The Green operator $\Gamma^0$ has an analytical expression in Fourier space, and hence the convolution $\Gamma^0 * \tau$ can be efficiently computed using FFT techniques [1, 2]. The $\Gamma^0$ operator thus defined enjoys a number of properties [7]:

$$ \Gamma^0 * \tau = 0 \iff \tau \in S \qquad (10) $$

$$ \langle \boldsymbol{e} : \mathbf{C}_0 : (\Gamma^0 * \boldsymbol{\tau}) \rangle_{\Omega} = \langle \boldsymbol{e} : \boldsymbol{\tau} \rangle_{\Omega} \text{ for all } \boldsymbol{\tau} \in \mathcal{L}_{\text{sym}}^2(\Omega) \text{ and } \boldsymbol{e} \in \mathbf{D}. \qquad (11) $$

To close this section, we recall the decomposition of any tensor field in $L_{\text{sym}}^2(\Omega)$ introduced in [17]; see also [19]. For any $t \in L_{\text{sym}}^2(\Omega)$, we define its $C_0$-norm as:

$$ \|t\|_{C_0} = \sqrt{\langle t : C_0 : t \rangle_{\Omega}}, \qquad (12) $$

and its decomposition:

$$ t^D = \Gamma^0 * (C_0 : t) \quad \text{and} \quad t^S = t - t^D. \qquad (13) $$

The following properties are easily established:

$$ t^D \in \mathbf{D}, \quad C_0 : t^S \in \mathbf{S} \quad \text{and} \quad \|t\|_{C_0}^2 = \|t^D\|_{C_0}^2 + \|t^S\|_{C_0}^2 \qquad (14) $$

and

$$ \langle t_1 : C_0 : t_2^S \rangle_{\Omega} = \langle t_1^S : C_0 : t_2 \rangle_{\Omega} \quad \text{and} \quad \langle t_1 : C_0 : t_2^D \rangle_{\Omega} = \langle t_1^D : C_0 : t_2 \rangle_{\Omega}, \qquad (15) $$

for all $t_1, t_2 \in L_{\text{sym}}^2(\Omega)$.

Hence, $t^D$ appears as the orthogonal projection of $t$ on $\mathbf{D}$ with respect to the $C_0$-norm. Unless $C_0$ is of the form $kI$, with $k$ a positive scalar and $I$ the fourth-order identity tensor operating on symmetric second order tensors ($I : t = t$ for any symmetric second order tensor $t$), the subspaces $\mathbf{D}$ and $\mathbf{S}$ are not orthogonal with respect to the $C_0$-norm, $t^S$ is not in $\mathbf{S}$ and it is not the orthogonal projection of $t$ on $\mathbf{S}$ with respect to the $C_0$-norm. $t^S$ is actually the projection (with respect to the $C_0$-norm) of $t$ on the subspace $S'$ of $t' \in L_{\text{sym}}^2(\Omega)$ such that $C_0 : t' \in \mathbf{S}$. Combining (10) with (13), we obtain the following characterizations of $\mathbf{S}$ and $\mathbf{D}$:

$$ (C_0^{-1} : \boldsymbol{\tau})^D = 0 \iff \boldsymbol{\tau} \in \mathbf{S}, \quad \text{and} \quad \boldsymbol{e}^S = 0 \iff \boldsymbol{e} \in \mathbf{D}. \qquad (16) $$

As a result, the $C_0$-norm of $(C_0^{-1} : \boldsymbol{\tau})^D$ can be considered as a measure of the equilibrium error for the stress field $\boldsymbol{\tau}$ and the $C_0$-norm of $\boldsymbol{e}^S$ can be considered as a measure of the compatibility error for the strain field $\boldsymbol{e}$.

Moreover, the residual $X$ can be decomposed as:

$$ X = X^S + X^D \text{ with } X^S = E - \boldsymbol{\epsilon}^S = (E - \boldsymbol{\epsilon})^S \text{ and } X^D = -\Gamma^0 * \boldsymbol{\sigma} = -(C_0^{-1} : \boldsymbol{\sigma})^D. \qquad (17) $$

Hence, $X = 0$ is equivalent to $\boldsymbol{\epsilon} - E \in \mathbf{D}$ and $\boldsymbol{\sigma} \in \mathbf{S}$.

It should be emphasized that multiplying $C_0$ by a strictly positive real does not change the above introduced decomposition (13).

### 3. AEM scheme for composites containing pores and rigid inclusions

Equation (6), which is the extension of the Lippmann-Schwinger to composites containing pores and rigid inclusions, is at the center of the so-called “FFT-based numerical homogenization techniques” first introduced by Moulinec and Suquet [1, 2]. The reader is referred to the recent review [16] for a more detailed description of this numerical technique.
The AEM scheme is started with $(\epsilon^0, \sigma^0) = (0,0)$ and its corresponding residual $X^0 = E$. For $n \ge 0$, if $X^n = 0$ then the pair $(\epsilon^n, \sigma^n)$ is a solution of the homogenization problem, and the iteration is stopped. Otherwise, compute $X^{n+1}$, $\epsilon^{n+1}$ and $\sigma^{n+1}$ as follows:

$$Z^n = \alpha : X^n + 2 ((I - \alpha) : X^n)^D \quad (18)$$

$$\lambda^n = \frac{\langle X^n : C_0 : Z^n \rangle_{\Omega}}{\langle Z^n : C_0 : Z^n \rangle_{\Omega}} \quad (19)$$

$$\epsilon^{n+1} = \epsilon^n + \lambda^n \alpha : X^n, \quad (20)$$

$$\sigma^{n+1} = \sigma^n + \lambda^n C_0 : (2I - \alpha) : X^n, \quad (21)$$

$$X^{n+1} = X^n - \lambda^n Z^n. \quad (22)$$

Here, the $\Omega$-periodic tensor field $\alpha(y)$ is defined over the whole unit cell as follows:

$$\forall y \in \Omega^r, \quad \alpha(y) = 0. \quad \forall y \in \Omega^p, \quad \alpha(y) = 2I. \quad \forall y \in \Omega^m, \quad \alpha(y) = 2(C(y) + C_0)^{-1} : C_0. \quad (23)$$

The polarization schemes correspond to the case where, instead of (19), the relaxation parameter $\lambda^n$ is set to a constant value, $1 - \gamma$, the damping parameter $\gamma$ being in the interval $[0, 1)$. The iteration of the Eyre-Milton scheme corresponds to $\gamma = 0$ ($\lambda^n = 1$) and the iteration of the augmented Lagrangian scheme corresponds to $\gamma = 0.5$ ($\lambda^n = 0.5$). In the AEM scheme, $\lambda^n$ is determined at each iteration $n$ such that as to minimize the $C_0$-norm of the residual at iteration $n+1$, $X^{n+1} = X^n - \lambda Z^n$, over all possible values of $\lambda$. Since the initialization in the original Eyre-Milton scheme is different from the zero-initialization $(\epsilon^0, \sigma^0) = (0,0)$, “$\gamma = 0$”-scheme will refer in the sequel to the AEM scheme (with zero-initialization) where the relaxation parameter $\lambda^n$ is set to 1 instead of (19).

Under a suitable geometric assumption on $\Omega^r$ and $\Omega^p$, and assuming that $C$ is uniformly coercive in the heterogeneous matrix, it has been proved in [17] that the sequence $(\epsilon^n, \sigma^n)$ linearly converges to the unique solution of the homogenization problem (4) such that: (a) the strain field in the pores is the unique solution to the following (Dirichlet) elasticity problem: the domain $\Omega^p$ is occupied by the reference material, there is no body forces and the prescribed displacements on its boundary are those of the heterogeneous matrix; (b) the stress field in the rigid inclusions is the unique solution to the following (Neumann) elasticity problem: the domain $\Omega^r$ is occupied by the reference material, there is no body forces and the prescribed tractions at its boundary are those of the heterogeneous matrix.

Numerical simulations have shown that the rate of convergence of the AEM scheme strongly depends on the choice of $C_0$, and it was recommended to choose $C_0$ of the same order as the matrix or the unknown homogenized material $C_{hom}$.

## 4. Accelerating AEM and “$\gamma = 0$” schemes

The idea is to update $C_0$ at each iteration with the twofold objective of speeding up convergence and making it independent of the choice of the initial value of $C_0$. We consider $C_0$ of the form $C_0 = k_0C_*$ where $C_*$ is a fixed given stiffness tensor and $k_0$ is a (strictly) positive multiplying constant which must be updated. The value of $k_0$ at iteration $n$ is noted $k_0^n$, and $\alpha^n$ denotes $\alpha$ when $C_0$ in (23) is set to $k_0^n C_*$.

The accelerated AEM scheme (A2EM) is started with $k_0^0 > 0$, $(\epsilon^0, \sigma^0) = (0,0)$ and $X^0 = E$. For $n \ge 0$, compute $\alpha^n$ with (23) and $X^{n+1}$, $\epsilon^{n+1}$, $\sigma^{n+1}$ and $k_0^{n+1}$ as follows:

$$Z^n = \alpha^n : X^n + 2 ((I - \alpha^n) : X^n)^D \quad (24)$$

$$\lambda^n = \frac{\langle X^n : C_* : Z^n \rangle_{\Omega}}{\langle Z^n : C_* : Z^n \rangle_{\Omega}} \quad (25)$$

$$\epsilon^{n+1} = \epsilon^n + \lambda^n \alpha^n : X^n, \quad (26)$$
$$ \boldsymbol{\sigma}^{n+1} = \boldsymbol{\sigma}^n + \lambda^n k_0^n \boldsymbol{C}_* : (2\boldsymbol{I} - \boldsymbol{\alpha}^n) : \boldsymbol{X}^n, \quad (27) $$

$$ \text{update } k_0^{n+1}, \quad (28) $$

$$ \boldsymbol{X}^{n+1} = \boldsymbol{X}^n - \lambda^n \boldsymbol{Z}^n + \left( (k_0^n)^{-1} - (k_0^{n+1})^{-1} \right) (\boldsymbol{C}_*^{-1} : \boldsymbol{\sigma}^{n+1})^D \quad (29) $$

We propose an updating rule for $k_0$ which results from the following analysis of the residual. Let $(\boldsymbol{\varepsilon}, \boldsymbol{\sigma}) \in \boldsymbol{L}_{\text{sym}}^2(\Omega) \times \boldsymbol{L}_{\text{sym}}^2(\Omega)$ be a pair of strain-stress fields which complies with the constitutive equations (3). It can be considered as a good approximation of the solution to the homogenization problem if its residual (7) is small enough. Thanks to the orthogonal decomposition of this residual (17), one can write:

$$ \|\boldsymbol{X}\|_{C_0}^2 = \|(E - \boldsymbol{\varepsilon})^S\|_{C_0}^2 + \|(\boldsymbol{C}_0^{-1} : \boldsymbol{\sigma})^D\|_{C_0}^2. \quad (30) $$

There are two contributions to the norm of the residual, $\|\boldsymbol{X}\|_{C_0}$: the error in compatibility, $\|(E - \boldsymbol{\varepsilon})^S\|_{C_0}$, and the error in equilibrium, $\|(\boldsymbol{C}_0^{-1} : \boldsymbol{\sigma})^D\|_{C_0}$. When these errors are not zero, their ratio depends on the choice of the reference medium $C_0 = k_0 C_*$: as $k_0$ increases, the relative contribution of the error in equilibrium to $\|\boldsymbol{X}\|_{C_0}$ decreases, and it increases as $k_0$ decreases. Since in the AEM scheme the strain is not compatible, nor the stress is in equilibrium, one has to ensure that both the projection of $E - \boldsymbol{\varepsilon} \neq 0$ on $S'$ and the projection of $\boldsymbol{C}_0^{-1} : \boldsymbol{\sigma} \neq 0$ on $D$ are small enough.

For this purpose, it is natural to normalize each projected field by its norm. We therefore define the normalized total error, $\Delta$, as:

$$ \Delta^2 = (\Delta^S)^2 + (\Delta^D)^2 \quad (31) $$

where

$$ 0 \le \Delta^S = \frac{\|(E - \boldsymbol{\varepsilon})^S\|_{C_0}}{\|E - \boldsymbol{\varepsilon}\|_{C_0}} = \frac{\|(E - \boldsymbol{\varepsilon})^S\|_{C_*}}{\|E - \boldsymbol{\varepsilon}\|_{C_*}} \le 1 \quad (32) $$

$$ 0 \le \Delta^D = \frac{\|(\boldsymbol{C}_0^{-1} : \boldsymbol{\sigma})^D\|_{C_0}}{\|\boldsymbol{C}_0^{-1} : \boldsymbol{\sigma}\|_{C_0}} = \frac{\|(\boldsymbol{C}_*^{-1} : \boldsymbol{\sigma})^D\|_{C_*}}{\|\boldsymbol{C}_*^{-1} : \boldsymbol{\sigma}\|_{C_*}} \le 1 \quad (33) $$

are the normalized compatibility and equilibrium errors, respectively.

Since $\Delta^S$ (respectively, $\Delta^D$) is a measure of the compatibility (respectively, equilibrium) error, then we set $\Delta^S = 0$ (respectively, $\Delta^D = 0$) for $\boldsymbol{\varepsilon} = E$ (respectively, $\boldsymbol{\sigma} = 0$). A natural stopping criterion of the AEM scheme (and its accelerated version) which ensures that both the projection of $E - \boldsymbol{\varepsilon} \neq 0$ on $S'$ and the projection of $\boldsymbol{C}_0^{-1} : \boldsymbol{\sigma} \neq 0$ on $D$ are small enough can be: $\Delta \le \text{tol}$ where $\text{tol}$ is a user-prescribed relative tolerance.

Remembering that the AEM scheme is formulated in terms of the residual $\boldsymbol{X}$, our idea is to choose $k_0$ in such a way that $\boldsymbol{X}$ gives equal weight to both normalized equilibrium and compatibility errors. We find that setting:

$$ k_0 = \frac{\|(\boldsymbol{C}_*^{-1} : \boldsymbol{\sigma})\|_{C_*}}{\|E - \boldsymbol{\varepsilon}\|_{C_*}}, \quad (34) $$

leads to the following equations:

$$ \Delta = \frac{\|\boldsymbol{X}\|_{C_0}}{\|E - \boldsymbol{\varepsilon}\|_{C_0}}, \quad \frac{\Delta^D}{\Delta} = \frac{\|\boldsymbol{X}^D\|_{C_0}}{\|\boldsymbol{X}\|_{C_0}} \quad \text{and} \quad \frac{\Delta^S}{\Delta} = \frac{\|\boldsymbol{X}^S\|_{C_0}}{\|\boldsymbol{X}\|_{C_0}}, \quad (35) $$

which means that the ratio of the compatibility and equilibrium errors in $\boldsymbol{X}$ is equal to the ratio of the corresponding normalized errors in $\Delta$.

Hence, we propose the following updating rule in (28), which we term *equiprojection rule*:

$$ k_0^{n+1} = \frac{\|(\boldsymbol{C}_*^{-1} : \boldsymbol{\sigma}^{n+1})\|_{C_*}}{\|E - \boldsymbol{\varepsilon}^{n+1}\|_{C_*}} \quad (36) $$
As pointed out in [20], the polarization-based schemes are actually related to a classical numerical solution method, the Douglas-Rachford splitting. Based on numerical simulations on heterogeneous linear and nonlinear materials (which may contain pores but no rigid inclusions), Schneider [21] recommended the use of an updating rule, initially proposed by Lorenz and Tran-Dinh [22] for general Douglas-Rachford splitting and Alternating Direction Method of Multipliers (ADMM) methods, given by:

$$C_* = I, \quad k_0^{n+1} = \frac{\|\boldsymbol{\sigma}^{n+1}\|_{L_{\text{sym}}^2}}{\|\boldsymbol{\varepsilon}^{n+1}\|_{L_{\text{sym}}^2}}. \qquad (37)$$

Finally, we will study also what we term the accelerated “$\gamma = 0$”-scheme which is A2EM scheme where the relaxation parameter $\lambda^n$ is set to 1 instead of (25).

## 5. Efficient implementation of the accelerated schemes

At first sight, A2EM scheme should require the use of an additional convolution at each iteration in order to update the residual (22). This is not the case if we organize the algorithm as follows. Let $(\boldsymbol{\varepsilon}, \boldsymbol{\sigma})$ be a pair of strain-stress fields complying to the constitutive equations (3) and $X$ its corresponding residual (7) with its decomposition (17). We have:

$$X^D = -(C_0^{-1} : \boldsymbol{\sigma})^D = -k_0^{-1} (C_*^{-1} : \boldsymbol{\sigma})^D = k_0^{-1} W \quad \text{where} \quad -(C_*^{-1} : \boldsymbol{\sigma})^D = W. \qquad (38)$$

The idea is to keep $W$ in memory and to update it without an additional convolution. Hence, the first modified A2EM is started with $k_0^0 > 0$, $(\boldsymbol{\varepsilon}^0, \boldsymbol{\sigma}^0) = (0, 0)$, $W^0 = 0$, $X^0 = E$ and $\Delta^0 = 1$. For $n \ge 0$, compute $X^{n+1}$, $W^{n+1}$, $\boldsymbol{\varepsilon}^{n+1}$, $\boldsymbol{\sigma}^{n+1}$, $k_0^{n+1}$ and $\Delta^{n+1}$ as follows:

$$Y^n = (\boldsymbol{\alpha}^n : X^n)^D \qquad (39)$$

$$Z^n = \boldsymbol{\alpha}^n : X^n - 2Y^n + 2(k_0^n)^{-1} W^n \qquad (40)$$

$$\lambda^n = \frac{\langle X^n : C_* : Z^n \rangle_{\Omega}}{\langle Z^n : C_* : Z^n \rangle_{\Omega}} \qquad (41)$$

$$\boldsymbol{\varepsilon}^{n+1} = \boldsymbol{\varepsilon}^n + \lambda^n \boldsymbol{\alpha}^n : X^n, \qquad (42)$$

$$\boldsymbol{\sigma}^{n+1} = \boldsymbol{\sigma}^n + \lambda^n k_0^n C_* : (2I - \boldsymbol{\alpha}^n) : X^n, \qquad (43)$$

$$W^{n+1} = (1 - 2\lambda^n) W^n + \lambda^n k_0^n Y^n, \qquad (44)$$

$$k_0^{n+1} = \frac{\|(C_*^{-1} : \boldsymbol{\sigma}^{n+1})\|_{C_*}}{\|E - \boldsymbol{\varepsilon}^{n+1}\|_{C_*}}, \qquad (45)$$

$$X^{n+1} = X^n - \lambda^n Z^n + ((k_0^n)^{-1} - (k_0^n)^{-1}) W^{n+1}, \qquad (46)$$

$$\Delta^{n+1} = \frac{\|X^{n+1}\|_{C_*}}{\|E - \boldsymbol{\varepsilon}^{n+1}\|_{C_*}}. \qquad (47)$$

There is only one convolution per iteration in the above scheme, namely for the computation of $Y$. As far as memory space is concerned, there are 5 variables to store and update in the above scheme: $X$, $W$, $\boldsymbol{\varepsilon}$, $\boldsymbol{\sigma}$ and the intermediate variable $Y$. Indeed, there is no need to store $Z$ which is only used in the computation of $\lambda$.

One can actually reduce to 4 variables to store if the auxiliary strain field $\boldsymbol{\varepsilon}_a$ introduced in [17] is used in the above scheme instead of both $\boldsymbol{\varepsilon}$ and $\boldsymbol{\sigma}$. Indeed, let $(\boldsymbol{\varepsilon}, \boldsymbol{\sigma})$ be a pair of strain-stress fields complying to the constitutive equations (3). Then, $\boldsymbol{\varepsilon}_a$ is given by:

$$\boldsymbol{\varepsilon}_a = \frac{1}{2} ((C_0)^{-1} : \boldsymbol{\sigma} + \boldsymbol{\varepsilon}), \qquad (48)$$
and one can easily check that $(\boldsymbol{\varepsilon}, \boldsymbol{\sigma})$ can be retrieved from $\boldsymbol{\varepsilon}_a$ by formulas:

$$ \boldsymbol{\varepsilon} = \boldsymbol{\alpha} : \boldsymbol{\varepsilon}_a, \quad \boldsymbol{\sigma} = \mathbf{C}_0 : (2\mathbf{I} - \boldsymbol{\alpha}) : \boldsymbol{\varepsilon}_a. \qquad (49) $$

Using these formulas at iteration $n$, we find that:

$$ \boldsymbol{\varepsilon}^{n+1} = \boldsymbol{\alpha}^n : (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n), \qquad (50) $$

$$ \boldsymbol{\sigma}^{n+1} = k_0^n \mathbf{C}_* : (2\mathbf{I} - \boldsymbol{\alpha}^n) (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n), \qquad (51) $$

$$ \boldsymbol{\varepsilon}_a^{n+1} = \frac{1}{2} \left( (k_0^{n+1})^{-1} (\mathbf{C}_*)^{-1} : \boldsymbol{\sigma}^{n+1} + \boldsymbol{\varepsilon}^{n+1} \right) = \frac{1}{2} \left( \frac{k_0^n}{k_0^{n+1}} (2\mathbf{I} - \boldsymbol{\alpha}^n) + \boldsymbol{\alpha}^n \right) : (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n). \qquad (52) $$

Hence, the second modified A2EM scheme is started with $k_0^0 > 0$, $\boldsymbol{\varepsilon}_a^0 = 0$, $\mathbf{W}^0 = 0$, $\mathbf{X}^0 = \mathbf{E}$ and $\Delta^0 = 1$. For $n \ge 0$, compute $\mathbf{X}^{n+1}$, $\mathbf{W}^{n+1}$, $\boldsymbol{\varepsilon}_a^{n+1}$, $k_0^{n+1}$ and $\Delta^{n+1}$ as follows:

$$ \mathbf{Y}^n = (\boldsymbol{\alpha}^n : \mathbf{X}^n)^D \qquad (53) $$

$$ \mathbf{Z}^n = \boldsymbol{\alpha}^n : \mathbf{X}^n - 2\mathbf{Y}^n + 2(k_0^n)^{-1} \mathbf{W}^n \qquad (54) $$

$$ \lambda^n = \frac{\langle \mathbf{X}^n : \mathbf{C}_* : \mathbf{Z}^n \rangle_{\Omega}}{\langle \mathbf{Z}^n : \mathbf{C}_* : \mathbf{Z}^n \rangle_{\Omega}} \qquad (55) $$

$$ k_0^{n+1} = k_0^n \frac{\|(2\mathbf{I} - \boldsymbol{\alpha}^n) : (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n)\|_{\mathbf{C}_*}}{\|\mathbf{E} - \boldsymbol{\alpha}^n : (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n)\|_{\mathbf{C}_*}}, \qquad (56) $$

$$ \boldsymbol{\varepsilon}_a^{n+1} = \frac{1}{2} \left( \frac{k_0^n}{k_0^{n+1}} (2\mathbf{I} - \boldsymbol{\alpha}^n) + \boldsymbol{\alpha}^n \right) : (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n), \qquad (57) $$

$$ \mathbf{W}^{n+1} = (1 - 2\lambda^n) \mathbf{W}^n + \lambda^n k_0^n \mathbf{Y}^n, \qquad (58) $$

$$ \mathbf{X}^{n+1} = \mathbf{X}^n - \lambda^n \mathbf{Z}^n + ((k_0^{n+1})^{-1} - (k_0^n)^{-1}) \mathbf{W}^{n+1}, \qquad (59) $$

$$ \Delta^{n+1} = \frac{\|\mathbf{X}^{n+1}\|_{\mathbf{C}_*}}{\|\mathbf{E} - \boldsymbol{\varepsilon}^{n+1}\|_{\mathbf{C}_*}}. \qquad (60) $$

The main steps are summarized in Algorithm 1.

<caption>**Algorithm 1** The A2EM algorithm with the equiprojection rule (36)</caption>

1: $\boldsymbol{\varepsilon} \leftarrow 0$
2: $\mathbf{X} \leftarrow \mathbf{E}$
3: $\mathbf{W} \leftarrow 0$
4: $k_0 \leftarrow k_0^0$
5: **repeat**
6: $\mathbf{Y} \leftarrow (\boldsymbol{\alpha} : \mathbf{X})^D$ where $\boldsymbol{\alpha}$ is computed with $\mathbf{C}_0 = k_0 \mathbf{C}_*$ and standard FFT are used with the analytical expression of the Green operator in Fourier space.
7: $\lambda \leftarrow \frac{\langle \mathbf{X} : \mathbf{C}_* : \mathbf{Z} \rangle_{\Omega}}{\langle \mathbf{Z} : \mathbf{C}_* : \mathbf{Z} \rangle_{\Omega}}$ with $\mathbf{Z} = \boldsymbol{\alpha} : \mathbf{X} - 2\mathbf{Y} + 2k_0^{-1}\mathbf{W}$
8: $\beta \leftarrow \|\mathbf{E} - \boldsymbol{\alpha} : (\boldsymbol{\varepsilon} + \lambda\mathbf{X})\|_{\mathbf{C}_*}$
9: $k_0^{\text{new}} \leftarrow k_0 \beta^{-1} \|(2\mathbf{I} - \boldsymbol{\alpha}) : (\boldsymbol{\varepsilon} + \lambda\mathbf{X})\|_{\mathbf{C}_*}$
10: $\boldsymbol{\varepsilon} \leftarrow \frac{1}{2} (\frac{k_0}{k_0^{\text{new}}} (2\mathbf{I} - \boldsymbol{\alpha}) + \boldsymbol{\alpha}) : (\boldsymbol{\varepsilon} + \lambda\mathbf{X})$
11: $\mathbf{X} \leftarrow \mathbf{X} - \lambda(\boldsymbol{\alpha} : \mathbf{X} - 2\mathbf{Y} + 2k_0^{-1}\mathbf{W}) + ((k_0^{\text{new}})^{-1} - (k_0)^{-1})((1 + 2\lambda)\mathbf{W} + \lambda k_0 \mathbf{Y})$
12: $\mathbf{W} \leftarrow (1 - 2\lambda)\mathbf{W} + \lambda k_0 \mathbf{Y}$
13: $k_0 \leftarrow k_0^{\text{new}}$
14: $\Delta \leftarrow \beta^{-1} \|\mathbf{X}\|_{\mathbf{C}_*}$
15: **until** Convergence criterion met $\Delta \le \text{tol}$
16: $\boldsymbol{\varepsilon} \leftarrow \boldsymbol{\alpha} : \boldsymbol{\varepsilon}$
17: **return** $\boldsymbol{\varepsilon}$
In the accelerated “$\gamma = 0$”-scheme, where $\lambda^n = 1$, the intermediate variables $\mathbf{Y}^n$ and $\mathbf{W}^n$ are not needed: there are only two variables to store: $\boldsymbol{\varepsilon}_a$ and $\mathbf{X}$. This scheme is started with $k_0^0 > 0$, $\boldsymbol{\varepsilon}_a^0 = 0$, $\mathbf{X}^0 = \mathbf{E}$ and $\Delta^0 = 1$. For $n \ge 0$, compute $\mathbf{X}^{n+1}, \boldsymbol{\varepsilon}_a^{n+1}, k_0^{n+1}$ and $\Delta^{n+1}$ as follows:

$$k_0^{n+1} = k_0^n \frac{\|(2\mathbf{I} - \boldsymbol{\alpha}^n) : (\boldsymbol{\varepsilon}_a^n + \lambda\mathbf{X}^n)\|_{C_*}}{\|\mathbf{E} - \boldsymbol{\alpha}^n : (\boldsymbol{\varepsilon}_a^n + \lambda\mathbf{X}^n)\|_{C_*}}, \quad (61)$$

$$\boldsymbol{\varepsilon}_a^{n+1} = \frac{1}{2} \left( \frac{k_0^n}{k_0^{n+1}} (2\mathbf{I} - \boldsymbol{\alpha}^n) + \boldsymbol{\alpha}^n \right) : (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n), \quad (62)$$

$$\mathbf{X}^{n+1} = \mathbf{E} - \boldsymbol{\alpha}^{n+1} : \boldsymbol{\varepsilon}_a^{n+1} - 2 ((\mathbf{I} - \boldsymbol{\alpha}^{n+1}) : \boldsymbol{\varepsilon}_a^{n+1})^D, \quad (63)$$

$$\Delta^{n+1} = \frac{\|\mathbf{X}^{n+1}\|_{C_*}}{\|\mathbf{E} - \boldsymbol{\alpha}^{n+1} : \boldsymbol{\varepsilon}_a^{n+1}\|_{C_*}}. \quad (64)$$

The main steps of the accelerated “$\gamma = 0$”-scheme are summarized in Algorithm 2.

<caption>**Algorithm 2** The accelerated “$\gamma = 0$”-scheme with the equiprojection rule (36)</caption>

1: $\boldsymbol{\varepsilon} \leftarrow 0$
2: $\mathbf{X} \leftarrow \mathbf{E}$
3: $k_0 \leftarrow k_0^0$
4: **repeat**
5: &nbsp;&nbsp;&nbsp;&nbsp; $\boldsymbol{\beta} \leftarrow \|\mathbf{E} - \boldsymbol{\alpha} : (\boldsymbol{\varepsilon} + \mathbf{X})\|_{C_*}$ where $\boldsymbol{\alpha}$ is computed with $C_0 = k_0 C_*$
6: &nbsp;&nbsp;&nbsp;&nbsp; $k_0^{\text{new}} \leftarrow k_0 \boldsymbol{\beta}^{-1} \|(2\mathbf{I} - \boldsymbol{\alpha}) : (\boldsymbol{\varepsilon} + \mathbf{X})\|_{C_*}$
7: &nbsp;&nbsp;&nbsp;&nbsp; $\boldsymbol{\varepsilon} \leftarrow \frac{1}{2} \left( \frac{k_0^{\text{new}}}{k_0^{\text{new}}} (2\mathbf{I} - \boldsymbol{\alpha}) + \boldsymbol{\alpha} \right) : (\boldsymbol{\varepsilon} + \mathbf{X})$
8: &nbsp;&nbsp;&nbsp;&nbsp; $k_0 \leftarrow k_0^{\text{new}}$
9: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{X} \leftarrow \mathbf{E} - \boldsymbol{\alpha} : \boldsymbol{\varepsilon} - 2((\mathbf{I} - \boldsymbol{\alpha}) : \boldsymbol{\varepsilon})^D$ where $\boldsymbol{\alpha}$ is computed with $C_0 = k_0 C_*$ and standard FFT are used with the analytical expression of the Green operator in Fourier space to compute the $D$-projection.
10: &nbsp;&nbsp;&nbsp;&nbsp; $\Delta \leftarrow \boldsymbol{\beta}^{-1} \|\mathbf{X}\|_{C_*}$
11: **until** Convergence criterion met $\Delta \le \text{tol}$
12: $\boldsymbol{\varepsilon} \leftarrow \boldsymbol{\alpha} : \boldsymbol{\varepsilon}$
13: **return** $\boldsymbol{\varepsilon}$

Both AEM and “$\gamma = 0$” schemes can be accelerated using the Lorenz and Tran-Dinh updating rule (37) instead of the equiprojection rule (36). In this case, equation (56) becomes:

$$k_0^{n+1} = k_0^n \frac{\|(2\mathbf{I} - \boldsymbol{\alpha}^n) : (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n)\|_{L_{\text{sym}}^2}}{\|\boldsymbol{\alpha}^n : (\boldsymbol{\varepsilon}_a^n + \lambda^n \mathbf{X}^n)\|_{L_{\text{sym}}^2}}, \quad (65)$$

and equation (61) becomes (65) where $\lambda^n$ is set to 1.

## 6. Illustrative applications

We test the performances of the AEM and “$\gamma = 0$” schemes and their accelerated versions, with the equiprojection updating rule (36) or the Lorenz and Tran-Dinh one (37), on doubly contrasted linear materials containing both pores and rigid inclusions. All the schemes have been implemented in the Janus Python library<fnref n="1" />. For the Green operator discretization scheme, we used the finite differences scheme proposed in [23] which results in improved accuracy of the local fields over more classical discretization strategies.

<footnote n="1">¹Available at https://github.com/sbrisard/janus.</footnote>
## 6.1. *The micro-structures*

We consider 2 micro-structures to compare the performances of the above mentioned schemes:

* *Micro-structure A:* a cubic unit cell consisting of a matrix (phase 1) with an infinitely rigid spherical inclusion of radius 0.25 (phase 2) located at its center and pores of radius 0.25 (phase 3) centered at the vertices of the unit cell (see Figure 1). The matrix is linear elastic and isotropic with Poisson ratio $v_1 = 0.3$ and shear modulus $\mu_1 = 1$. The unit cell is discretized with a 32 × 32 × 32 grid.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** A. 3-phase 3D medium with double contrast, of size 32*32*32. The pore is in light blue, the inclusion in grey and the matrix is hidden.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** B. Matrix-fiber-pore system, of size 256*256*256 (the matrix is hidden).</figcaption></figure>

* *Micro-structure B:* a fiber composite with gas injection. The gas is injected in the fiber-reinforced composites to reduce the weight and some effects like residual stress. This results in a three-phase micro-structure containing identical spherical pores (5% volume fraction) and fibers with random orientation (10% volumetric fraction). See Figure 2. The fibers are infinitely rigid and the stiffness in the pores is null. In order to avoid any
artificial contact, there is a minimal distance of $4\sqrt{3}$ voxels (corresponding to 4 voxels in diagonal) between two fibers, two pores or between a fiber and a pore. The unit cell is discretized with a $256 \times 256 \times 256$ grid.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** $\Delta$ at iterations 20(a), 50(b), 100(c) and 500(d) as a function of $k_0$. Reference material is of the form $C_0 = k_0C_{matrix}$. Comparison between AEM scheme and “$\gamma = 0$” scheme on micro-structure A.</figcaption></figure>

## 6.2. Influence of the reference material on the convergence of the AEM and “$\gamma = 0$”-schemes

We consider microstructure A with reference material of the form $C_0 = k_0C_{matrix}$, which means that we set $C_* = C_{matrix}$. Figure 3 represents the normalized total error $\Delta$ at iterations 20, 50, 100 and 500 as a function of $k_0$. It clearly shows that the rate of convergence is optimal for some value of $k_0$, $k_0^{opt}$, close to 1, as predicted by the estimation given in [17], and that the rate of convergence is getting much worse as $k_0$ deviates from $k_0^{opt}$. Another important observation is that the rate of convergence of both schemes (AEM and “$\gamma = 0$”) is the same for $k_0$ smaller than $k_0^{opt}$, and that AEM converges faster than “$\gamma = 0$” otherwise.

## 6.3. Influence of the initial reference material on the convergence of the accelerated schemes

We still consider microstructure A with reference material of the form $C_0 = k_0^n C_{matrix}$ where $k_0^n$ is given by the updating rule (36) or by the Lorenz and Tran-Dinh one (37). Figure 4 (respectively,
Figure 5) represents the normalized total error $\Delta$ at iterations 20, 50, 100 and 500 of the A2EM scheme (respectively, accelerated “$\gamma = 0$”-scheme) as a function of the initial value of $k_0$, $k_0^0$, for both updating rules. In addition, we have also plotted the AEM scheme (respectively, the “$\gamma = 0$”-scheme) for which $k_0$ remains constant (equal to $k_0^0$) for all iterations.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** $\Delta$ at iterations 20(a), 50(b), 100(c) and 500(d) as a function of $k_0^0$, the initial value of $k_0$. The reference medium is of the form $C_0 = k_0^n C_{matrix}$. Comparison on micro-structure A between the AEM scheme ($k_0^n = k_0^0$ remains unchanged) and the A2EM scheme with the equiprojection updating rule (36) and the Lorenz and Tran-Dinh updating rule (37).</figcaption></figure>

Figures 4 and 5 show that $\Delta$ almost does not depend on $k_0^0$ after some iterations (50 in this case) for both accelerated schemes with both updating rules. More precisely, if $k_0^0$ is less than $k_0^{opt}$ (the optimal value of $k_0$ for the non accelerated scheme) then the performance of the accelerated scheme is very comparable to non accelerated schemes with $k_0 = k_0^{opt}$. On the other hand, for $k_0^0$ greater than $k_0^{opt}$, the rate of convergence of the accelerated scheme is a bit slower than non accelerated schemes with $k_0 = k_0^{opt}$, but much faster than non accelerated schemes with $k_0 = k_0^0$.

After some iterations (here between 20 and 50), $k_0^n$ reaches an asymptotic value which is independent of the initial value $k_0^0$, as shown in Figure 6. Moreover, it is seen that: (a) both updating rules lead to comparable rates of convergence, and (b) the accelerated AEM converges faster than the accelerated “$\gamma = 0$” as illustrated in Figure 7 which represents $\Delta$ at iterations 20, 50, 100 and 500 as a function of $k_0^0$ for both AEM and “$\gamma = 0$” schemes accelerated with the updating rule (36).
Finally, we present in Figure 8 the normalized total error $\Delta$ at iterations 50, 100, 500 and 1000 of the AEM and “$\gamma = 0$”-schemes and their accelerated versions with the equiprojection updating rule (36) as a function of $k_0^0$, the initial value of $k_0$ (which is not updated for AEM and “$\gamma = 0$”-schemes). It is seen that the convergence rate for micro-structure B is slower than the convergence rate for the simpler micro-structure A, probably because the stress, strain fields are less smooth. Again AEM is faster than “$\gamma = 0$” and there is a common optimal value of reference material (around 10 times the matrix stiffness and 6 times the homogenized material) for which the convergence is remarkably fast. This optimal performance is retrieved by A2EM and accelerated “$\gamma = 0$” schemes for any choice of $k_0^0$, and especially for $k_0^0$ lower than the optimal value.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** $\Delta$ at iterations 20(a), 50(b), 100(c) and 500(d) as a function of $k_0^0$, the initial value of $k_0$. The reference medium is of the form $C_0 = k_0^n C_{\text{matrix}}$. Comparison on micro-structure A between the “$\gamma = 0$”-scheme ($k_0^n = k_0^0$ remains unchanged) and the accelerated “$\gamma = 0$”-scheme with the equiprojection updating rule (36) and the Lorenz and Tran-Dinh updating rule (37).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** $k_0^n$ as a function of the initial value $k_0^0$, for $n = 20$(a) and $n = 50$(b), in the A2EM scheme for both updating rules. Micro-structure A.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** $\Delta$ at iterations 20(a), 50(b), 100(c) and 500(d) as a function of $k_0^0$, the initial value of $k_0$. The reference medium is of the form $C_0 = k_0^n C_{matrix}$. Comparison on micro-structure A between the A2EM and accelerated "$\gamma = 0$" schemes with the equiprojection updating rule (36).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** $\Delta$ at iterations 50(a), 100(b), 500(c) and 1000(d) as a function of $k_0^0$, the initial value of $k_0$. The reference medium is of the form $C_0 = k_0^n C_{matrix}$. Comparison on micro-structure B between the AEM and "$\gamma = 0$" schemes and their accelerated version with the equiprojection updating rule (36).</figcaption></figure>

## 7. Conclusions

We proposed to accelerate the AEM scheme, a recent iterative method for the numerical solution of the Lippmann-Schwinger equation with periodic boundary conditions for the case of linear materials containing pores and rigid inclusions. This scheme is derived from the classical Eyre-Milton scheme [6] and the polarization-based schemes [8], where at each iteration, the *direction* of the increment is preserved, while its *amplitude* is optimized. When the amplitude is not optimized, the scheme is called "$\gamma = 0$"-scheme. The main contributions of the present paper are summarized below.

(1) Numerical simulations with AEM and "$\gamma = 0$" schemes on micro structures containing both pores and rigid inclusions have been conducted with reference material of the form $C_0 = k_0 C_*$, where $C_*$ is a fixed given stiffness tensor and $k_0$ is a (strictly) positive multiplying constant. It is observed that the convergence of these schemes is very fast for some optimal value of $k_0$, $k_0^{opt}$. Besides, the rate of convergence is getting much worse as $k_0$ deviates from $k_0^{opt}$. Our simulations show that the optimal reference material could deviate from the estimations given in [17]: the optimal reference material is 10 times the matrix and 6 times the homogenized material for micro-structure B.
(2) The rate of convergence of AEM and “$γ = 0$” are the same for $k_0$ lower than $k_0^{opt}$ while AEM is faster than “$γ = 0$” otherwise.

(3) We proposed to accelerate AEM and “$γ = 0$” schemes for $k_0 \neq k_0^{opt}$ by updating the reference material of the form $C_0 = k_0^n C_{matrix}$ where $k_0^n$ at iteration $n$ is given by the proposed equiprojection updating rule (36) or by the Lorenz and Tran-Dinh one (37).

(4) Efficient implementation of these accelerated schemes in terms of memory and number of convolutions per iteration has been proposed. It is shown that only one convolution per iteration is needed for both accelerated schemes, and that A2EM needs twice memory as accelerated “$γ = 0$”.

(5) Numerical simulations with both A2EM and accelerated “$γ = 0$” schemes and both updating rules on micro structures containing pores and rigid inclusions show that, after some iterations (50 in the considered example), the rate of convergence almost does not depend on the choice of the initial value of $k_0$ and that it is comparable to the rate of convergence of the non accelerated schemes with the optimal value $k_0 = k_0^{opt}$.

(6) The numerical simulations show that the A2EM scheme is faster than the accelerated “$γ = 0$” scheme.

In summary, we recommend the use of A2EM and accelerated “$γ = 0$” schemes with the equiprojection rule, and to start these schemes with a reference material that underestimates the homogenized one. In future works, these schemes will be applied to nonlinear materials, where their robustness with respect to the material contrast will be particularly helpful. There is also a need for additional mathematical studies on the convergence of the proposed accelerated schemes.

**CRediT authorship contribution statement**

* Martin Dolbeau: data curation, software, validation, visualization.

* Jérémy Bleyer: data curation, software, validation, visualization.

* Karam Sab: conceptualization, formal analysis, investigation, methodology, writing.

**Data availability**

The scripts used to produce the results presented in the paper can be found at [24, 10.5281/zenodo.12706218].

**Declaration of interests**

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

**References**

[1] H. Moulinec and P. Suquet, “A Fast Numerical Method for Computing the Linear and Nonlinear Mechanical Properties of Composites”, C. R. Acad. Sci., Sér. IIA Earth Planet. Sci. **318** (1994), no. 11, pp. 1417–1423.

[2] H. Moulinec and P. Suquet, “A Numerical Method for Computing the Overall Response of Nonlinear Composites with Complex Microstructure”, Comput. Methods Appl. Mech. Eng. **157** (1998), no. 1-2, pp. 69–94.

[3] J. Korringa, “Theory of Elastic Constants of Heterogeneous Media”, J. Math. Phys. **14** (1973), no. 4, pp. 509–513.

[4] R. Zeller and P. H. Dederichs, “Elastic Constants of Polycrystals”, Phys. Status Solidi B Basic Res. **55** (1973), no. 2, pp. 831–842.
[5] E. Kröner, "On the Physics and Mathematics of Self-Stresses", in *Topics in Applied Continuum Mechanics* (J. L. Zeman and F. Ziegler, eds.), Springer: Vienna, 1974, pp. 22–38.

[6] D. J. Eyre and G. W. Milton, "A Fast Numerical Scheme for Computing the Response of Composites Using Grid Refinement", *Eur. Phys. J. AP* **6** (1999), no. 01, pp. 41–47.

[7] J. C. Michel, H. Moulinec and P. Suquet, "A Computational Scheme for Linear and Non-Linear Composites with Arbitrary Phase Contrast", *Int. J. Numer. Methods Eng.* **52** (2001), no. 1-2, pp. 139–160.

[8] V. Monchiet and G. Bonnet, "A Polarization-Based FFT Iterative Scheme for Computing the Effective Properties of Elastic Composites with Arbitrary Contrast", *Int. J. Numer. Methods Eng.* **89** (2012), no. 11, pp. 1419–1436.

[9] H. Moulinec and F. Silva, "Comparison of Three Accelerated FFT-based Schemes for Computing the Mechanical Response of Composite Materials", *Int. J. Numer. Methods Eng.* **97** (2014), no. 13, pp. 960–985.

[10] J. Zeman, J. Vondřejc, J. Novák and I. Marek, "Accelerating a FFT-based Solver for Numerical Homogenization of Periodic Media by Conjugate Gradients", *J. Comput. Phys.* **229** (2010), no. 21, pp. 8065–8071.

[11] S. Brisard and L. Dormieux, "FFT-based Methods for the Mechanics of Composites: A General Variational Framework", *Comput. Mater. Sci.* **49** (2010), no. 3, pp. 663–671.

[12] L. Gélébart and R. Mondon-Cancel, "Non-Linear Extension of FFT-based Methods Accelerated by Conjugate Gradients to Evaluate the Mechanical Behavior of Composite Materials", *Comput. Mater. Sci.* **77** (2013), pp. 430–439.

[13] M. Kabel, T. Böhlke and M. Schneider, "Efficient Fixed Point and Newton–Krylov Solvers for FFT-based Homogenization of Elasticity at Large Deformations", *Comput. Mech.* **54** (2014), no. 6, pp. 1497–1514.

[14] M. Schneider, "An FFT-based Fast Gradient Method for Elastic and Inelastic Unit Cell Homogenization Problems", *Comput. Methods Appl. Mech. Eng.* **315** (2017), pp. 846–866.

[15] M. Schneider, "On the Barzilai–Borwein Basic Scheme in FFT-based Computational Homogenization", *Int. J. Numer. Methods Eng.* **118** (2019), no. 8, pp. 482–494.

[16] M. Schneider, "A Review of Nonlinear FFT-based Computational Homogenization Methods", *Acta Mech.* **232** (2021), no. 6, pp. 2051–2100.

[17] K. Sab, J. Bleyer, S. Brisard and M. Dolbeau, "An FFT-based adaptive polarization method for infinitely contrasted media with guaranteed convergence", *Comput. Methods Appl. Mech. Eng.* **427** (2024), article no. 117012.

[18] G. W. Milton, *The Theory of Composites*, Cambridge University Press, 2002.

[19] C. Bellis and P. Suquet, "Geometric variational principles for computational homogenization", *J. Elasticity* **137** (2019), pp. 119–149.

[20] M. Schneider, D. Wicht and T. Böhlke, "On Polarization-Based Schemes for the FFT-based Computational Homogenization of Inelastic Materials", *Comput. Mech.* **64** (2019), no. 4, pp. 1073–1095.

[21] M. Schneider, "On non-stationary polarization methods in FFT-based computational micromechanics", *Int. J. Numer. Methods Eng.* **122** (2021), no. 22, pp. 6800–6821.

[22] D. A. Lorenz and Q. Tran-Dinh, "Non-stationary Douglas–Rachford and alternating direction method of multipliers: adaptive step-sizes and convergence", *Comput. Optim. Appl.* **74** (2019), no. 1, pp. 67–92.

[23] F. Willot, "Fourier-based schemes for computing the mechanical response of composites with accurate local fields", *C. R. Méc. Acad. Sci. Paris* **343** (2015), no. 3, pp. 232–245.

[24] M. Dolbeau, K. Sab and J. Bleyer, *Supplementary material to 'Accelerating the Adaptive Eyre–Milton FFT-based method for infinitely contrasted media'*, Zenodo, 2024.