Isabelle Ramière
email: isabelle.ramiere@cea.fr

Thomas Helfer
email: thomas.helfer@cea.fr

Iterative residual-based vector methods to accelerate fixed point iterations

Keywords: Fixed point method, Acceleration iterative technique, Vector sequences, Residual-based method, Nonlinear mechanics, Mechanical behaviours

Fixed point iterations are still the most common approach to dealing with a variety of numerical problems such as coupled problems (multi-physics, domain decomposition,...) or nonlinear problems (electronic structure, heat transfer, nonlinear mechanics, ...). Methods to accelerate fixed point iteration convergence or more generally sequence convergence have been extensively studied since the 1960's. For scalar sequences, the most popular and efficient acceleration method remains the ∆ 2 of Aitken. Various vector acceleration algorithms are available in the literature, which often aim at being multi-dimensional generalizations of the ∆ 2 method.

In this paper, we propose and analyze a generic residual-based formulation for accelerating vector sequences. The question of the dynamic use of this residualbased transformation during the fixed point iterations for obtaining a new accelerated fixed point method is then raised. We show that two main classes of such iterative algorithms can be derived and that this approach is generic in that various existing acceleration algorithms for vector sequences are thereby recovered.

In order to illustrate the interest of such algorithms, we apply them in the field of nonlinear mechanics on a simplified "point-wise" solver used to perform mechanical behaviour unit testings. The proposed test cases clearly demonstrate that

Introduction

Fixed point iterations (also called Picart iterations) are used in many applications to deal with nonlinear equations. In fact, it is well-known that any nonlinear system of equations F (X) = 0 with F : R N → R N is equivalent to a fixed point problem G(X) = X with G : R N → R N . The fixed point iterations method is therefore the simplest way to obtain nonlinear solutions without any a priori knowledge about F , its derivative in particular is not required. On the other hand, fixed point iterations converge generally only linearly and very slowly. This is the reason why acceleration algorithms are required. As the fixed point iterations method generates a sequence, classical sequence acceleration methods can be applied. Many of such acceleration or extrapolation algorithms for vector sequences were proposed in the 20 th century [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF], especially in the 1960's (e.g. [START_REF] Shanks | Non linear transformations of divergent and slowly convergent sequences[END_REF][START_REF] Wynn | Acceleration techniques for iterated vector and matrix problems[END_REF][START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF][START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF]) and in the 1990's (e.g. [START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF][START_REF] Graves-Morris | Extrapolation method for vector sequences[END_REF][START_REF] Sedogbo | Some convergence acceleration processes for a class of vector sequences[END_REF][START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF]). Most of them are generalizations of the powerful Aitken's ∆ 2 scalar acceleration algorithm [START_REF] Aitken | On Bernouilli's numerical solution of algebraic equations[END_REF].

Applications of such vector acceleration techniques for the fixed point iterations method are various and diverse. For example, fixed point acceleration methods have recently enjoyed a renewed interest in the development of advanced numerical methods (domain decomposition [START_REF] Marini | A relaxation procedure for domain decomposition methods using finite elements[END_REF][START_REF] Garbey | Acceleration of the Schwarz method for elliptic problems[END_REF], multigrid [START_REF] Washio | Krylov subspace acceleration for nonlinear multigrid schemes[END_REF][START_REF] Jemcov | Acceleration and Stabilization of Algebraic Multigrid Solver Applied to Incompressible Flow Problems[END_REF],...) and of multiphysics coupling in a black-box context (fluid-structure interaction [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Degroote | Vierendeels, Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction[END_REF], thermomechanics [START_REF] Erbts | Accelerated staggered coupling schemes for problems of thermoelasticity at finite strains[END_REF][START_REF] Michel | Simulation of pellet-cladding interaction with the PLEIADES fuel performance software environment[END_REF], ...).

The aim of this article is to propose a common framework to build new vector acceleration methods. This formalism is build as a generalization of the nonlinear hybrid acceleration procedure [START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF] but can also be viewed as an extension of both the reduced rank extrapolation method [START_REF] Eddy | Extrapolation to the limit of a vector sequence[END_REF] and the Anderson's method [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF]. By the way, some of the most popular and efficient extrapolation methods can be recovered using this generic approach.

As no function derivative is involved in the proposed approach, this unified formalism can be applied on various type of problems as for example black box solver systems. Moreover, by construction, high-dimensionality problems can also be treated within this formalism.

The paper is organized as follows. We first propose in Section 2 a survey of most popular vector acceleration methods in a common fixed point sequence framework. In particular, we present them as iterative methods that consist in dynamically generating a new accelerated fixed point sequence with the transformed iterates. Then, we introduce in Section 3 an extension of the nonlinear hybrid acceleration procedure [START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF] in order to deal with linear combinations of fixed point iteration residuals. We show in Section 4 that two main classes of iterative (or dynamic) vector acceleration methods can be derived from this generic residual-based formulation. Various already known acceleration methods are retrieved within these two classes. Finally, Section 5 is dedicated to some performance comparisons between acceleration methods stem from these two formulations. The application concerns a simplified "point-wise" solver, developed within the PLEIADES platform [START_REF] Plancq | PLEIADES : a unified environment for multidimensional fuel performance modeling[END_REF] as part of the MFront software [START_REF] Cea | MFront Web Site[END_REF], which is mainly used to perform mechanical behaviour unit testings. This solver allows us, among other things, to compare our accelerated algorithms to the standard Newton-Raphson algorithm used by default in most mechanical finite element solvers.

Survey of iterative vector acceleration methods

The aim of this section is to briefly describe most of the existing vector sequence acceleration methods. As it is pointed out in the literature (see for example [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF][START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF]), static sequence transformations from an existing sequence are less efficient than dynamic or cycling sequence transformations that consist in considering the accelerated iterate as the new iterate of the sequence. Therefore, we shall only present here the iterative version of the vector acceleration methods.

Moreover, for the sake of clarity, we shall focus on fixed point iteration sequences as they are the most common sequences dealt with in an engineering context (even the Newton's method generates a fixed point sequence). Within this fixed point iteration framework, expressions of existing vector sequence acceleration approaches can be more easily compared.

As a universal acceleration method for all types of sequences cannot exist (see for example [START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF]), in this paper we shall only consider transformations that have been designed to accelerate first order (or linearly converging) sequences. Only few authors proposed some approaches to accelerate sequences of order greater than one, the interested reader is referred to [START_REF] Jones | Extrapolation for higher orders of convergence[END_REF][START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF].

Scalar sequences

Let us begin this survey be the scalar case, which is the basis of most of the vector acceleration methods. Considering a scalar nonlinear fixed point equation

x = g(x), x ∈ R, g : R → R, (1)
the most popular and powerful acceleration method of the basic fixed point substitution iteration

x n+1 = g(x n ) (2)
remains the ∆ 2 of Aitken [START_REF] Aitken | On Bernouilli's numerical solution of algebraic equations[END_REF] and its recursive application by the Steffensen algorithm which lead to a second-order method (e.g. [START_REF] Henrici | Elements of Numerical Analysis[END_REF][START_REF] Nievergelt | Aitken's and Steffensen's accelerations in several variables[END_REF]):

x n+1 = x n - (g(x n ) -x n ) 2 g(g(x n )) -2g(x n ) + x n = x n - (∆x n ) 2 ∆ 2 x n (3)
with ∆ the difference operator,

∆x n = g(x n ) -x n , ∆g(x n ) = g(g(x n )) -g(x n ) and ∆ 2 x n = ∆g(x n ) -∆x n .
Various equivalent formulations exist, including for example

x n+1 = g(x n ) - ∆x n ∆g(x n ) ∆ 2 x n (4)
or

x n+1 = g(g(x n )) - (∆g(x n )) 2 ∆ 2 x n . ( 5
)
Whatever the expression, each Steffensen iteration requires two new basic fixed point iterations (or function evaluations) to be applied. The secant method (see [START_REF] Henrici | Elements of Numerical Analysis[END_REF] for example) can also be viewed as a ∆ 2 transformation requiring only one new fixed point iteration

x n+1 = x n - (x n -x n-1 )(g(x n ) -x n ) (g(x n ) -x n ) -(g(x n-1 ) -x n-1 ) = x n - (x n -x n-1 )∆x n ∆x n -∆x n-1 (6)
or equivalently

x n+1 = g(x n ) - (g(x n ) -g(x n-1 ))∆x n ∆x n -∆x n-1 . (7)
The secant method reduces to the Steffensen method by applying the acceleration only at every other step of the basic fixed point iteration.

The secant method converges with the order 1+ √ 5 2

(see for example [START_REF] Henrici | Elements of Numerical Analysis[END_REF]), but as it requires only one new basic fixed point iterate (or one new function evaluation per iteration), its efficiency index defined as r 1/s with r the order of the method and s the number of function evaluations per iteration required by the method, is little better

( 1+ √ 5 2 1.618) than the Steffensen's one ( √ 2 1.414).
Another famous acceleration algorithm is the relaxation method (see [START_REF] Quarteroni | Numerical Mathematics[END_REF] for example), sometimes called mixing algorithm [START_REF] Eyert | A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences[END_REF][START_REF] Fang | Two classes of multisecant methods for nonlinear acceleration[END_REF]:

x n+1 = ωg(x n ) + (1 -ω)x n . ( 8
)
The determination of the optimal relaxation parameter ω opt is the main drawback of this algorithm because it depends on the derivative of g which is a priori unknown. Hence empiric static values are often set in practice (most often ω ≡ 0.5, see for example [START_REF] Erbts | A partitioned solution approach for electro-thermo-mechanical problems[END_REF]). Using the equivalent formulation

x n+1 = x n + ω(g(x n ) -x n ) (9)
it is obvious that the secant and the Steffensen methods are particular relaxation methods with a varying relaxation parameter. Hence, in the literature the secant method is also called Aitken relaxation or dynamic relaxation [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Erbts | Accelerated staggered coupling schemes for problems of thermoelasticity at finite strains[END_REF].

Many other scalar acceleration techniques, also called extrapolation techniques, have been proposed since Aitken, but these methods often aim to be generalizations of the Aitken's process, see [START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF] for a detailed description of some of these algorithms. In particular, let us mention the -algorithm [START_REF] Wynn | On a device for computing the e m (S n ) transformation[END_REF] n

-1 = 0, n 0 = x n , n+1 0 = g(x n ) n k+1 = n+1 k-1 + 1 n+1 k -n k , (10)
whose 2 version, i.e. x n+1 = n 2 , is equivalent to the Steffensen method. An interesting generalization of several acceleration methods is proposed by the E-algorithm [START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF].

Vector sequences

Practical problems deal most often with many unknowns and systems of nonlinear equations. So let us now consider the vector fixed point equation

X = G(X), X ∈ R N , G : R N → R N , N ∈ N , N > 1, (11)
leading to the basic vector fixed point sequence

X n+1 = G(X n ). (12)
Vector sequence acceleration methods have been extensively studied in the literature [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF][START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF][START_REF] Wynn | Acceleration techniques for iterated vector and matrix problems[END_REF][START_REF] Smith | Extrapolation methods for vector sequences[END_REF][START_REF] Graves-Morris | Extrapolation method for vector sequences[END_REF][START_REF] Sedogbo | Some convergence acceleration processes for a class of vector sequences[END_REF][START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF][START_REF] Fang | Two classes of multisecant methods for nonlinear acceleration[END_REF]. Most of them are extensions of scalar acceleration methods with the following definition of the inverse of a real vector

X -1 = X X 2 . ( 13
)
The norm involved in the preceding definition is often chosen to be the Euclidean norm even if other choices can be made, see [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF] for example. More generally, the inverse of a vector can be defined via the inner product (denoted here by • ) with another vector

X -1 = Y Y • X , ∀Y ∈ R N . ( 14
)
As the ∆ 2 -Steffensen method can be expressed in different equivalent manners in the scalar case (see for example Eqs (3)-( 5)), many vector formulations based on the definition ( 13) and ( 14) are available. For example the method proposed by Irons and Tuck in [START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF] is obtained from Eq. ( 5) with the definition (13) of the inverse of a vector:

X n+1 = G(G(X n )) - ∆G(X n ) • ∆ 2 X n ∆ 2 X n 2 ∆G(X n ) (15)
where as previously

∆X n = G(X n ) -X n , ∆G(X n ) = G(G(X n )) -G(X n ) and ∆ 2 X n = ∆G(X n ) -∆X n .
Sometimes in the literature, the following expression is also attributed to Irons and Tuck (see for example [START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF][START_REF] Sedogbo | Some convergence acceleration processes for a class of vector sequences[END_REF]) but it seems in fact to be the approach of Lemaréchal [31]

X n+1 = X n - ∆X n • ∆ 2 X n ∆ 2 X n 2 ∆X n . (16)
This method can be directly obtained from the combination of definition [START_REF] Washio | Krylov subspace acceleration for nonlinear multigrid schemes[END_REF] and the Steffensen's scalar equation [START_REF] Wynn | Acceleration techniques for iterated vector and matrix problems[END_REF].

The following methods are based on the second definition of the inverse of a vector [START_REF] Jemcov | Acceleration and Stabilization of Algebraic Multigrid Solver Applied to Incompressible Flow Problems[END_REF]. Starting from the scalar equation ( 4), the approach of Graves-Morris [START_REF] Graves-Morris | Extrapolation method for vector sequences[END_REF] gives

X n+1 = G(X n ) - ∆X n 2 ∆X n • ∆ 2 X n ∆G(X n ), (17)
while the first step of the vector A-algorithm of Sedogbo [START_REF] Sedogbo | Some convergence acceleration processes for a class of vector sequences[END_REF] can be expressed as

X n+1 = G(X n ) - ∆G(X n ) 2 ∆G(X n ) • ∆ 2 X n ∆X n . (18)
Starting from the scalar Equation ( 5), the method of Jennings [START_REF] Jennings | Accelerating the convergence of matrix iterative processes[END_REF] writes

X n+1 = G(G(X n )) - ∆X n • ∆G(X n ) ∆X n • ∆ 2 X n ∆G(X n ), (19)
as Zienkiewicz and Lohner [START_REF] Zienkiewicz | Accelerated relaxation or direct solution. Future prospects for FEM[END_REF] proposed the following vector acceleration method

X n+1 = G(G(X n )) - ∆G(X n ) 2 ∆G(X n ) • ∆ 2 X n ∆G(X n ). (20)
The vector -algorithm of Wynn [START_REF] Wynn | Acceleration techniques for iterated vector and matrix problems[END_REF] is also a direct extension of the scalaralgorithm (see [START_REF] Wynn | On a device for computing the e m (S n ) transformation[END_REF] and Eq. ( 10)) with the definition (13) of the inverse of a vector. Hence the vector 2 -algorithm writes

X n+1 = G(X n ) + ∆G(X n ) ∆X n 2 -∆X n ∆G(X n ) 2 ∆ 2 X n 2 . ( 21
)
In his article, Macleod [START_REF] Macleod | Acceleration of vector sequences by multi-dimensional ∆ 2 methods[END_REF] compared 9 multi-dimensional ∆ 2 methods that reduce in one dimension to the Steffensen method. Some of them are listed previously but others were also introduced. In particular, the application of the scalar Steffensen algorithm on each component of the vector was tested.

However the results presented in [START_REF] Macleod | Acceleration of vector sequences by multi-dimensional ∆ 2 methods[END_REF], in accordance with [START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF], show that vector extrapolation algorithms (based on projection approaches) are more interesting than the scalar ones on each component. In these test cases, the Irons and Tuck [START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF] approach was the more efficient.

Many variants of the multi-dimensional secant method [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Fang | Two classes of multisecant methods for nonlinear acceleration[END_REF] can be formally obtained in the same manner from the scalar equations ( 6) and [START_REF] Graves-Morris | Extrapolation method for vector sequences[END_REF].

Other vector extensions of scalar extrapolation methods are detailed in [START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF][START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF].

These methods can be viewed as multi-step generalizations of iterative vector ∆ 2 approaches, in the sense that they require more basic fixed point iterations.

The nonlinear hybrid approach proposed in [START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF] introduces the vector ∆ k method:

X n+1 = X n - ∆X n • ∆ k+1 X n ∆ k+1 X n 2 ∆ k X n (22)
where

                         ∆ k X n = ∆ k-1 G(X n ) -∆ k-1 (X n ) ∆ k G l (X n ) = ∆ k-1 G l+1 (X n ) -∆ k-1 G l (X n ) ∆ 0 = Id G l (X n ) = G(G(...G l times (X n ))) G 0 = Id .
For k = 1, the Lemaréchal's method is recovered, see Eq. ( 16) . This general-ization for k > 1 seems suitable to deal with bifurcation problems (see [START_REF] Marder | A bifurcation problem in E-layer equilibria[END_REF] for example).

Another interesting generalization of multi-dimensional secant approaches has been originally proposed by Anderson [4]. This method has been directly built for vector sequences and is based, like the nonlinear hybrid approach [START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF], on a minimization process. This approach is a one-step method as it is applied every fixed point iteration. As it has been pointed out in some recent articles [START_REF] Eyert | A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences[END_REF][START_REF] Fang | Two classes of multisecant methods for nonlinear acceleration[END_REF][START_REF] Walker | Anderson acceleration for fixed-point iterations[END_REF], this method seems to be efficient in many situations. In the following sections, we shall see that the reduced rank extrapolation method of Eddy [START_REF] Eddy | Extrapolation to the limit of a vector sequence[END_REF] can be viewed as a multi-step Anderson's method.

A generic residual-based acceleration approach

Let (X n ) n be a sequence of vectors of R N converging to an unknown limit X, M an integer such that M ≤ N and (Z i n ) n , i = 1, .., M arbitrary vector sequences of R N converging to zero. We define the transformation Y that transforms the sequence

(X n ) n into a new sequence (Y n ) n Y n = X n - M i=1 λ i n Z i n ( 23
)
with λ i n , i = 1, .., M , some scalars.

When M = 1, expression (23) reduces to the nonlinear hybrid procedure of Brezinski and Chehab [START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF]. Obviously if the λ i n are constant parameters independent of n : λ i n ≡ λ i , ∀n, the sequence (Y n ) n converges to X and hence the transformation Y is regular. As in practice the λ i n vary with n, the regularity of Y can not be ensured for all convergent sequences. This drawback is quite widespread in the existing extrapolation methods, it is the case for example of the Aitken's ∆ 2 process (e.g. [START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF]).

However, under the usual assumption (see for example [START_REF] Aitken | On Bernouilli's numerical solution of algebraic equations[END_REF][START_REF] Brezinski | Extrapolation Methods: Theory and Practice[END_REF][START_REF] Jbilou | Vector extrapolation methods. Applications and numerical comparison[END_REF]) that consists in setting the same coefficients in the expressions of Y n and Y n+1 , an estimation of X can be obtained by minimizing δY n = Y n+1 -Y n with respect to these coefficients. In our case, we can write

Y n+1 = X n+1 - M i=1 λ i n Z i n+1 , (24)
and hence

δY n = δX n - M i=1 λ i n δZ i n ( 25
)
where the forward difference operator δ acts on the lower index:

δS n = S n+1 - S n .
By construction, the optimal parameters λ i n , i = 1, .., M lead to

190 δY n ≤ δX n (26)
for the chosen norm. The vector sequence (Y n ) n therefore converges faster than

(X n ) n .
For the Euclidean norm, an expression of the λ i n can be obtained using the normal equation solution 1 . Defining the following N × M matrices

Z n = (Z 1 n . . . Z M n ), δZ n = (δZ 1 n . . . δZ M n )
the least-square minimization gives

λ n =      λ 1 n . . . λ M n      = (δZ T n δZ n ) -1 δZ T n δX n . (27)
In practice, if a vector δZ i n is collinear to another vector δZ j n , j = i, the matrix 195 (δZ T n δZ n ) becomes singular and then its inversion is impossible. In this case, the corresponding λ i n is set to zero and the system is reduced.

The expression ( 27) of λ n is then used to define two sequence transformations which will lead to extrapolation methods (see Section 4):

Y n = X n -Z n (δZ T n δZ n ) -1 δZ T n δX n (28)
or

Y n+1 = X n+1 -Z n+1 (δZ T n δZ n ) -1 δZ T n δX n . (29)
The generic formulation [START_REF] Erbts | A partitioned solution approach for electro-thermo-mechanical problems[END_REF] in the end is very close to the polynomial extrapolation method formalism as described in [START_REF] Jbilou | Vector extrapolation methods. Applications and numerical comparison[END_REF]. The reduced rank extrapolation method [START_REF] Eddy | Extrapolation to the limit of a vector sequence[END_REF] is recovered for Z i n = δX n+i-1 which reduces for M = 1 to the Lemaréchal's method [START_REF]Une méthode de résolution de certains systèmes non linéaires bien posés[END_REF] (see Eq. ( 16)).

The choice Z i n = δX n-i enables us to deal with a linear combination of the previous sequences residuals. In this case, expression [START_REF] Erbts | A partitioned solution approach for electro-thermo-mechanical problems[END_REF] leads for M = 1 to a vector extension of the Steffensen's method of Eq. ( 4) whereas a generalization of the Irons and Tuck method (see [START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF] or Eq. ( 15)) is obtained with expression [START_REF] Wynn | On a device for computing the e m (S n ) transformation[END_REF]. Moreover, the basic idea of Anderson [4] is also recovered with this approach. For this choice, the kernel of the transformation Y is the same as the Shanks [START_REF] Shanks | Non linear transformations of divergent and slowly convergent sequences[END_REF] transformation:

k i=0 a i (X n-i -X) = 0 (30)
with k i=0 a i = 0 and a 0 .a k = 0.

Many other residuals Z i n could be chosen, for example [START_REF] Brezinski | Nonlinear hybrid procedures and fixed point iterations[END_REF] for the ∆ k method (see Eq. ( 22)).

Z i n = δ k X n-i , k > 1 as proposed in

Two classes of iterative residual-based acceleration methods

This section is devoted to the application of the residual-based transformations described in Section 3 to accelerate the fixed point iteration method (see Eq. ( 12)). As mentioned in the introduction of Section 2, we are only interested in iterative acceleration methods. Moreover, we want to obtain one-step iterative methods (i.e. only one new fixed point evaluation between each new accelerate iterate) as opposed to multi-step iterative methods such as the Irons and Tuck method [START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF] (see Eq. ( 15)), the vector -algorithm [START_REF] Wynn | Acceleration techniques for iterated vector and matrix problems[END_REF] (see Eq. ( 21)) or the polynomial extrapolation algorithms performed in [START_REF] Jbilou | Vector extrapolation methods. Applications and numerical comparison[END_REF]. Considering that the iterative acceleration process applied to the fixed point iteration generates two sequences (X n ) n and (G(X n )) n (see [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF] for example), two main formalisms of iterative methods are available. They differ in the way the two coupled sequences are taken into account in the definition of Y n and Y n+1 in Equation ( 23)- [START_REF] Nievergelt | Aitken's and Steffensen's accelerations in several variables[END_REF].

To illustrate our point, for each formalism we shall decline the iterative method obtained for the most common choice Z i n = δX n-i in the transformation (29) of Section 3. Methods derived from other choices of Z i n can be easily deduced from this example.

As it is not guaranteed that the transformations Y defined in Section 3 are regular with respect to the fixed point limit X, the convergence criterion must still focus on the standard fixed point iteration residual ∆X n = G(X n ) -X n to obtain the desired solution.

First class : crossed sequences method

The first formalism focuses on the basic fixed point sequence

(G(X n )) n Y n = G(X n-1 ) - M i=1 λ i n Z i n (31) Y n+1 = G(X n ) - M i=1 λ i n Z i n+1 , (32)
and takes into account the accelerated sequence (X n ) n in the definition of the vanishing sequences

(Z i n ) n , i = 1, ..M . Indeed, Z i n depends in this case on the fixed point iteration residual r(X n ) ∆X n = G(X n ) -X n . ( 33
)
For example, the case Z i n = δX n-i in the generic approach of Section 3 becomes for the iterative crossed sequences method

Z i n = ∆X n-i as Z i n = δ 2 X n-i-1 be- comes Z i n = ∆(X n-i ) -∆(X n-i-1 ) = G(X n-i ) -X n-i -G(X n-i-1 ) + X n-i-1 .
Example For the case Z i n = ∆X n-i , the crossed sequences method gives

δY n = (G(X n ) -G(X n-1 )) - M i=1 λ i n (∆X n-i+1 -∆X n-i ). (34)
Using the notation

R n = (∆X n-1 • • • • • • ∆X n-M ) δR n = ((∆X n -∆X n-1 ) • • • • • • (∆X n+1-M -∆X n-M )),
we obtain

X n+1 = Y n+1 = G(X n ) -R n+1 (δR T n δR n ) -1 δR T n (G(X n ) -G(X n-1 )). (35)
For M = 1, the standard vector secant acceleration or dynamic relaxation method is recovered (e.g. [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Erbts | Accelerated staggered coupling schemes for problems of thermoelasticity at finite strains[END_REF]). This method corresponds formally to the vector extension of the scalar secant method [START_REF] Graves-Morris | Extrapolation method for vector sequences[END_REF] with the definition (13) of the inverse of the vector:

X n+1 = G(X n ) - (G(X n ) -G(X n-1 )) • (∆X n -∆X n-1 ) ∆X n -∆X n-1 2 ∆X n . (36)
In a manner similar to the scalar case, if this acceleration method is applied alternately with a basic fixed point iteration step, the Irons and Tuck method [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF] is recovered.

Second class: alternate sequences method

In this case, Y n is concerned by the sequence

(X n ) n as Y n+1 is concerned by the sequence (G(X n )) n : Y n = X n - M i=1 λ i n Z i n ( 37
) Y n+1 = G(X n ) - M i=1 λ i n Z i n+1 , (38)
where Z i n is a function of X n-i , i = 0, ..M whereas Z i n+1 depends only on G(X n-i ), i = 0, ..M . We may have for example

   Z i n = δX n-i = X n-i+1 -X n-i Z i n+1 = δG(X n-i ) = G(X n-i+1 ) -G(X n-i ) or    Z i n = δ 2 X n-i-1 = X n-i+1 -2X n-i + X n-i-1 Z i n+1 = δ 2 G(X n-i-1 ) = G(X n-i+1 ) -2G(X n-i ) + G(X n-i-1 )
.

For any choice of

Z i n = δ k X n-i-k+1 , k > 1, the expression of δY n = Y n+1 -Y n is then a linear combination of fixed point residuals ∆X n = G(X n ) -X n , the
coefficients of which are identical to the sequence coefficients in the accelerated sequences equations [START_REF] Jbilou | Vector extrapolation methods. Applications and numerical comparison[END_REF] or [START_REF] Verpeaux | Algorithmes et méthodes du code éléments finis Cast3M[END_REF]. In this case, if the acceleration process is chosen to be a two-step method (the acceleration is then applied every other step of the basic fixed point iteration), the equation δY n to minimize becomes the same in both formalisms. However, except for the case k = 1 and M = 1, the expression of X n+1 remains different between each formalism, and so the obtained acceleration methods differ.

Remark When a nonlinear equation F (X) = 0 is under consideration, the residual F (X n ) is often used instead of G(X n )-X n , see for example [START_REF] Eyert | A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences[END_REF][START_REF] Fang | Two classes of multisecant methods for nonlinear acceleration[END_REF][START_REF] Verpeaux | Algorithmes et méthodes du code éléments finis Cast3M[END_REF].

Example For the case Z i n = δX n-i / Z i n+1 = δG(X n-i ), the alternate sequences method reads

δY n = (G(X n ) -X n ) - M i=1 λ i n δ(G(X n-i ) -X n-i ) = ∆X n - M i=1 λ i n δ∆X n-i = ∆X n - M i=1 λ i n (∆X n-i+1 -∆X n-i ). (39)
With the same definition of δR n as in the crossed sequence method and introducing

δG n = (δG(X n-i ) • • • δG(X n-M )),
the accelerate iterate writes:

X n+1 = Y n+1 = G(X n ) -δG n (δR T n δR n ) -1 δR T n ∆X n . (40)
This extrapolation method was first introduced by Anderson in the 1960's [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF] and has recently attracted more attention [START_REF] Eyert | A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences[END_REF][START_REF] Fang | Two classes of multisecant methods for nonlinear acceleration[END_REF][START_REF] Walker | Anderson acceleration for fixed-point iterations[END_REF] 2 . Moreover, this approach is similar to the so-called "interface quasi-Newton method" [START_REF] Degroote | Vierendeels, Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction[END_REF] or "reduced order models" method used in [START_REF] Vierendeels | Implicit coupling of partitioned flui-structure interaction problems with reduced order models[END_REF].

For M = 1, this method reduces to another vector secant method

X n+1 = G(X n ) - (∆X n -∆X n-1 ) • ∆X n ∆X n -∆X n-1 2 (G(X n ) -G(X n-1 )). (41)
Both vector secant methods of Eq. ( 36) and Eq. ( 41) can be obtained directly from Eq. ( 7) with the definition of the inverse of a vector (13) but they differ from the numerator's term involved in the inner product. Moreover, these two secant methods reduce to the Irons and Tuck method if the acceleration method is applied at every other step.

Application

MTest open-source tool

To illustrate the performances of some of the residual-based acceleration algorithms presented in this paper (see Section 4), we introduced them in an open-source tool named MTest which is developed in the PLEIADES plateform [START_REF] Plancq | PLEIADES : a unified environment for multidimensional fuel performance modeling[END_REF] as part of the MFront project [START_REF] Cea | MFront Web Site[END_REF]. MFront is a code generator that enables engineers and researchers in the field of structural mechanics to implement various types of mechanical behaviour easily, reliably and efficiently.

MTest has been designed to test mechanical behaviours on a single material point by imposing constraints on each component of the strains or the stresses 2 In the original formulation, Anderson proposed to define X n+1 as a relaxation between Y n+1 and Yn. However, in practice the relaxation parameter was set to 1 and hence X n+1 ≡ Y n+1 .

σ. Since and σ are rank-2 symmetric tensors, they will be described by two vectors E and Σ. In 3D the number of components of E and Σ is equal to

N = 6.
For the sake of clarity, in the sequel we shall only describe the case in which all the stress components are imposed. Imposing strain components is accomplished by introducing Lagrange multipliers which increase the problem size.

Equilibrium equation. At each time step between t and t+∆ t, a strain increment X = E| t+∆ t -E| t which respects the mechanical equilibrium is sought. This strain increment thus satisfies the following nonlinear equation :

R (X) = Σ| t+∆ t (X, ∆ t, V | t ) -Σ| imp t+∆ t = 0, (42)
where Σ| t+∆ t is a generally nonlinear function computing the stresses as a result of the mechanical behaviour integration, V is a set of internal state variables known at the beginning of the time step and whose value at the end of the time step is also a result of the mechanical behaviour integration, and Σ| imp t+∆ t are the imposed stress values. In this paper, the mechanical behaviour integration will be considered as a black box.

On the use of the Newton-Raphson algorithm. Problem (42) can be solved by the Newton-Raphson algorithm if the behaviour provides the so-called consistent tangent operator ∂ Σ| t+∆ t ∂X [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF]. This strategy is considered to be the default in most implicit finite element solvers dedicated to structural mechanics [START_REF]Abaqus Analysis User's Manual[END_REF][START_REF] Ansys | USER material subroutine USERMAT[END_REF][START_REF] Edf | Algorithm for nonlinear quasi-static analysis (operator STAT_NON_LINE)[END_REF][START_REF]Northwest Numerics and Modeling[END_REF]. The Cast3M finite element solver is a noticeable exception and its default algorithm paved the way to the method presented in this section [START_REF]CEA, Cast3M Web Site[END_REF][START_REF] Verpeaux | Algorithmes et méthodes du code éléments finis Cast3M[END_REF].

Providing consistent tangent operator is only feasible analytically in certain cases when the behaviour integration is done using an implicit scheme [START_REF] Besson | Numerical implementation of constitutive models[END_REF].

Even in this case, its computation may require a tremendous amount of work.

It is also worth noting that quadratic convergence of the Newton-Raphson algorithm strongly depends on the quality of this operator and can easily be lost due to implementation mistakes or, even worse, to a poor convergence criteria for the behaviour integration. This point is specifically outlined in the Abaqus user guide [START_REF]Abaqus Analysis User's Manual[END_REF].

In structural mechanics, the use of the consistent tangent operator may have other drawbacks: this matrix may be unsymmetric and not definite positive (due to softening for example), which forbids the use of standard efficient linear solvers.

Alternative algorithm. It is therefore worth considering algorithms not relying on the consistent tangent operator to solve Equation [START_REF] Ansys | USER material subroutine USERMAT[END_REF]. One interesting characteristic of mechanical behaviour is that it usually introduces an elastic operator K which allows us to rewrite Equation ( 42) as a fixed point problem using a quasi-Newton method:

X = G (X) with G (X) = X -K -1 R (X) . (43)
In practice, for the great majority of mechanical behaviours, the elastic operator K is "stiff enough" to guarantee that the spectral radius of the Jacobian matrix of G is less than 1 and hence that the fixed point iterations method will converge.

Solving Equation ( 43) with a standard fixed point iterations method thus has the following advantages:

-K is easy to compute.

-K is symmetrical and definite positive.

-K can be factorized only once through the whole computation.

In practice, if finite strains are considered or if the elastic properties evolve with external parameters, such as temperature, it may be worth updating this matrix from time to time.

However, as mentioned in the introduction of this paper, such a fixed point iterations method leads to a very slow linear convergence. Therefore, use of acceleration procedures is mandatory to obtain an efficient method.

Stopping criterion and accuracy. Whatever the algorithm used, iterations are stopped once the absolute maximum norm of the residual The results exhibit two linear stages. At the very beginning, the behaviour is elastic. The initial slope has a magnitude close to the Young modulus, which is equal to 145.2 GPa (quasi vertical line on Figure 1). At the end of the test, the stress is also almost a linear function of the strain: the effective slope is around 300 smaller than the initial one. This ratio gives an estimation of how close the elastic operator is from the consistent tangent operator. Performances of the standard fixed point iterations for solving the quasi-Newton problem [START_REF] Edf | Algorithm for nonlinear quasi-static analysis (operator STAT_NON_LINE)[END_REF] are expected to decrease as this ratio increases.

∆ X n = G(X n ) -X n
Furthermore, the behaviour is integrated in this case with an implicit scheme that supplies the consistent tangent operator: we can thus compare the accelerated fixed point algorithms proposed in this paper to the standard Newton-Raphson algorithm.

Polycrystals test case

The behaviour used in this second test case is the result of the Berveiller-Zaoui homogenization scheme to a FCC-polycristal made of 30 grains [START_REF] Berveiller | An extension of the self-consistent scheme to plastically-flowing polycrystals[END_REF]. The resulting behaviour is orthotropic. FCC crystals have 12 sliding systems, the plasticity in each grain is described by flow rules, no by sliding systems, derived from dislocation dynamics (see the Code-Aster documentation for details [START_REF] Edf | Mono and polycrystalline elastoviscoplastic constitutive laws[END_REF]).

This behaviour introduces 1272 state variables, which is very large and unusual.

The material is submitted to a uniaxial tensile test: the axial strain ε zz monotonically and linearly increases from 0 at the beginning of the test to 5 % at the end of the test. The loading is described in 15 time steps. All stress components are null, except the axial stress σ zz . The result of this test case is reported in Figure 2. As for the Chaboche test case, Figure 2 shows two linear stages. The slope ratio between those stages in this case is close to 70. The standard fixed point iterations method is hence expected to converge faster than in the Chaboche test case.

Implicit schemes are not considered for the integration of this behaviour for two main reasons:

1. the Jacobian matrix of the implicit system would be fairly complex to compute;

2. such an implicit scheme would require to solving a 1272 × 1272 dense system at each Newton-Raphson iteration.

The integration is thus performed using a standard Runge-Kutta algorithm. As a consequence, the consistent tangent operator is not available, so the equilibrium Equation ( 42) cannot be solved using a Newton-Raphson algorithm.

The solver must then rely on standard fixed point iterations based on Equation [START_REF] Edf | Algorithm for nonlinear quasi-static analysis (operator STAT_NON_LINE)[END_REF]. Acceleration methods for fixed point iterations are thus extremely appealing for this test case.

Results

Acceleration methods compared

In this section we shall compare different residual-based acceleration algorithms generated thanks to the two formalisms detailed in Section 4. As the dimension of the involved problem is N = 6 (see Section 5.1), we will compare acceleration methods involving only two or three iterations.

Concerning the methods involving the last two iterations, we study the performances of the residual methods with M = 1 and Z 1 n = δX n-1 :

• the crossed secant method or crossed 1-δ method (often called Aitken relaxation or dynamic relaxation in the literature [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF])

X n+1 = G(X n ) - (G(X n ) -G(X n-1 )) • (∆X n -∆X n-1 ) ∆X n -∆X n-1 2 ∆X n ; (44)
• the alternate secant method or alternate 1-δ method (also called Anderson extrapolation method with M = 1 [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF])

X n+1 = G(X n ) - (∆X n -∆X n-1 ) • ∆X n ∆X n -∆X n-1 2 (G(X n ) -G(X n-1 )); (45)
• the well-known Irons and Tuck method [START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF] in which both preceding methods reduce if applied alternately with a standard fixed point iteration

X n = G(X n-1 ) X n+1 = G(X n ) - ∆X n • ∆ 2 X n-1 ∆ 2 X n-1 2 ∆X n . (46)
However, this is no longer true for the crossed 2-δ method. So we can also consider the following variant of the crossed 2-δ method

• crossed 2-δ bis method X n+1 = G(X n ) -λ 1 n (G(X n ) -X n ) -λ 2 n (G(X n ) -X n-1 ), with λ 1 n , λ 2 n minimizing δY n = (G(X n ) -G(X n-1 )) -λ 1 n (∆X n -∆X n-1 ) -λ 2 n (G(X n ) -G(X n-1 ) -X n-1 + X n-2 ). (51)
In the following sections, we shall compare the acceleration performances of each residual methods in terms of number of iterations required for the whole loading (total number of iterations) or at each time step to reach the convergence criterion. We designate by number of iterations the number of evaluations of the fixed point function G. As the behaviour integration (used to evaluate R, see Section 5.1) is the most costly part of the computation per iteration, this measure is relevant. Others situations are discussed in Section 5.3.4.

Two iteration residual methods comparison

Chaboche test case. For the Chaboche test case, Table 1 presents the total number of iterations (sum of all the time steps) required by the two iteration residual methods introduced in Section 5.3.1 as well as by the standard fixed point iterations method. As the consistent tangent operator is available for this test case, we also report in the table the total number of iterations required by the standard second-order Newton method. Moreover, the numerical order of convergence obtained for each method is also given.

The first conclusion to be drawn from Table 1 is that all the listed acceleration methods converge to the fixed point solution in this case. We can also note the very poor convergence of the fixed point iterations method. case, the alternate secant method seems the most powerful two iteration residual method to accelerate the fixed point convergence whatever the precision required. Moreover, the numbers of iterations involved by the alternate secant method are quite close to those of the Newton's method, which confirms the efficiency and the performance of this method. Numerically, the convergence order of both secant methods is around 1.6 which is in good agreement with the theoretical order of convergence of the scalar secant approach (see Section 2.1).

However the rate of convergence in this case is in favour of the alternate secant method. For the Irons and Tuck approach, no numerical order of convergence appears clearly during the simulations.

Figure 3 presents a detailed view of the number of iterations required per acceleration method and per time step for a desired precision of ε = 10 -8 .

This figure confirms that whatever the time step and hence the loading, the alternate secant method is the more efficient two iteration residual-based method.

The behaviour of this acceleration method is furthermore very stable and similar to the Newton's one. From this figure, we can also remark that the Irons and Tuck method has some convergence difficulties in the final linear stage. If we take a deep look at time step #10 for example, we can see in Figure 4 that when the Irons and Tuck method is applied, the residual decreases slowly after a precision of 10 -7 (but still faster than for the fixed point iteration method).

This figure also shows the quite chaotic residual convergence of the crossed se- cant method and the Newton's method. For the Newton's method, we get back to the situation where the initial solution is far from the converged solution and hence the convergence of the method isn't monotone. In this case, the Newton's method may diverge and is then often combined with line-search methods for example.

Polycrystals test case. Table 2 recapitulates the total number of iterations required by every two iteration residual-based method for the Polycrystals test case. The total number of iterations of the standard fixed point iteration method is also reported. This test case is not suitable to the evaluation of the numerical order of convergence: in fact no method presents a regular order of convergence, for the three approaches, it ranges from 0.5 to 3 according to the time step.

As for the Chaboche test case, all the listed acceleration methods converge to For this test case, the best acceleration method depends on the required precision. On the whole, the crossed secant method seems to be the most efficient, especially for the most stringent precision. The Irons and Tuck approach is also interesting in this case. On Figure 5, the number of iterations required by each acceleration method is represented versus the time step for ε = 10 -8

(Figure 5 Conclusion. These two test cases enable us to conclude that both vector secant methods are efficient in accelerating the fixed point iterations method convergence. According to the Chaboche test case, they also exhibit the same order of convergence of 1.6, which is in good agreement with the theoretical scalar secant convergence order. Finally, the best acceleration factor (or rate of convergence) between both approaches seems to depend on the problem under consideration.

It appears that the more slowly the fixed point iteration method converges, the more the alternate secant method is efficient compared to the crossed secant method.

Three iteration residual methods comparison

In this section, we present the performances obtained for the two test cases with the five 'three iteration residual methods' listed in section 5.3.1, Eqs. ( 47)- [START_REF] Edf | Mono and polycrystalline elastoviscoplastic constitutive laws[END_REF]. For each method and for each time step, at iteration #2, the corresponding secant method is first applied.

Chaboche test case. The total numbers of iterations required for the Chaboche test case by the three iteration methods studied in this paper are reported in The first conclusion to be drawn from this table is that all acceleration approaches except the crossed 2-δ bis method, converge to the fixed point solution. However, as shown in Figure 7, the accelerated sequence generated by the crossed 2-δ bis method converges, but not towards the limit of the fixed point iteration method. This example clearly demonstrates the fact that in order to guarantee that the solution obtained thanks to an acceleration algorithm is truly the desired solution of the fixed point problem, the convergence criterion must remain based on the fixed point residual ∆X n = G(X n ) -X n whatever the acceleration algorithm performed.

The second conclusion is that the crossed δ 2 method is not a good acceleration method, since the number of iterations and the order of convergence obtained with this method are similar to those of the fixed point iterations method, see As for the Chaboche test case, the crossed 2-δ bis method does not converge to 580 the fixed point solution, even if the accelerated sequence converges. For this test case, the crossed δ 2 approach is worse in term of number of iterations than the standard fixed point iteration method (see Table 2). Moreover for ε = 10 -12 , this method has not yet reached the convergence criterion at time step #1 after 10.000 iterations!

As opposed to the Chaboche test case, here the crossed 2-δ method is on the whole little more efficient than the alternate δ 2 method but less interesting than the crossed secant method (see Table 2). Here again, it is the alternate 2-δ method that proves to be the most reliable acceleration method involving the last three iterations. This method is in particular more efficient than the crossed secant, which was the most powerful two iterations method for this test case. These conclusions are clearly summarized by Figure 8.

Conclusion.

For the test cases under consideration here, the acceleration methods based on the δ 2 residual do not seem relevant whatever the formalism. 600 Secondly, for both test cases, the alternate 2-δ, already known as the Anderson extrapolation method with M = 2 in literature, is the best outlined iterative acceleration method. It reduces at least to the most efficient vector secant method and even may improve it. This last conclusion is in good agreement with the experiments of Anderson [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF] in which the alternate δ methods were successfully 605 applied on a special class of nonlinear equations.

Discussion

In this paragraph we shall discuss some points related to the application of the proposed acceleration methods to problems with a large number N of degrees of freedom.

It is worth noting that the cost specifically attached to the acceleration methods proposed is generally negligible in comparison to the cost of the evaluation of the fixed point function. For example, in the case of the crossed secant method (44), we must compute two scalar products and update a vector. The cost of those operations grows linearly with N and can be easily parallelized by highly scalable implementations (such as BLAS routines, for example). This remark still holds true in distributed memory.

Moreover, the conclusion that the Newton-Raphson method, when available, appears as the best nonlinear resolution method is difficult to generalize as N grows. Indeed, in this case, due to the cost of Jacobian matrix factorization (which is negligible in our test cases), the number of iterations is no longer a reliable indicator of the efficiency of this method.

Concerning the memory space, the conclusions drawn from previous numerical examples encourage us to consider alternate M-δ methods with M being greater than 2. However it is worth underlining that in this case, M + 1 iterates must be stored and this may limit the performances of the approach for high dimensional problems. In those cases, a balance would have to be found between acceleration factor and memory space.

Conclusions and perspectives

In this paper, we have proposed and analyzed a generic residual-based acceleration approach to construct sequence acceleration processes. This approach can be viewed as a generalization of various existing vector acceleration methods (also called extrapolation methods). Then, this generic formalism has been derived in two main classes of methods to accelerate fixed point iterations convergence: the crossed and the alternate sequences approaches.

We proposed to test the performances of the obtained acceleration algorithms on the first order fixed point iterations method obtained from the simulation of mechanical behaviour unit testings when the consistent tangent coherent oper-ator is unknown or difficult to obtain. A set of iterative residual methods based on the two or three last iterates has been compared. It follows that some of the proposed acceleration methods, especially those derived from the alternate approach, are really efficient methods: they are able for example to compete against the second-order Newton-Raphson method.

In the future, it would be worthwhile to confirm the performances of such strategies on structural mechanics solvers involving large numbers of degrees of freedom.

Finally, as the proposed iterative residual-based approach is really generic and independent of the application context, the proposed formalism and the derived algorithms can be useful for many other applications.

σFigure 1 :

 1 Figure 1: Chaboche test case

Figure 2 :

 2 Figure 2: Polycrystals_dd_cfc test case

Figure 3 :

 3 Figure 3: Number of iterations per time stepε = 10 -8 -Chaboche test case -Two iteration residual methods

Figure 4 :

 4 Figure 4: Residuals for time step #10ε = 10 -8 -Chaboche test case -Two iteration residual methods

  (a)) and ε = 10 -12 (Figure 5(b)).

Figure 5 :

 5 Figure 5: Number of iterations per time step -Polycrystals test case -Two iteration residual methods

1 Figure 6 :

 16 Figure 6: Residuals for the Polycrystals test case -Two iteration residual methods

Figure 7 :

 7 Figure 7: Residuals of the crossed 2-δ bis method for ε = 10 -8 for time step #6 -Chaboche test case

Figure 8 :

 8 Figure 8: Number of iterations per time step for ε = 10 -8 -Polycrystals test case -Three iteration residual methods

Figure 9

 9 Figure 9 confirms that crossed extrapolation methods have a more chaotic convergence behaviour per time step than the alternate methods and the fixed point iteration method. Moreover, the efficiency of the alternate 2-δ method as well as

Figure 9 :

 9 Figure 9: Residuals for the time step #5 and ε = 10 -8 -Polycrystals test case -Three iteration residual methods

Table 1 :

 1 Total number of iterations and numerical convergence order for the Chaboche test case -Two iteration residual methods

	The residual-

Table 2 :

 2 Total number of iterations for the Polycrystals test case -Two iteration residual methods the problem solution whatever the precision required. Moreover, the standard fixed point iteration method converges definitely faster than in the Chaboche case. This is in good agreement with the elastic operator properties in this case, see Section 5.2.2. However, the proposed residual-based acceleration methods still accelerate the fixed point iterations convergence.

	510

Table 3 .

 3 The numerical orders of convergence are also reported in this table.

		crossed δ 2 alternate δ 2 crossed 2-δ alternate 2-δ crossed 2-δ bis
	ε=10 -6	9509	83	80	70	71
	ε=10 -8	19098	102	145	85	Not Cvg
	ε=10 -12	37357	122	687	102	Not Cvg
	Cvg order	1	-	1.6	1.6	-

Table 3 :

 3 Total number of iterations and numerical convergence order for the Chaboche test case -Three iteration residual methods

Table 1 .

 1 The other three iteration methods (alternate δ 2 method, crossed 2-δ and alternate 2-δ) are interesting, the alternate methods being the most efficient, especially more powerful than the crossed secant method. No regular order of convergence had been numerically obtained for the alternate δ 2 method.

	It varies from 0.5 to 2.5 according to the time step. It is worth noting that the
	2-δ methods recover the order of convergence of the corresponding secant (or
	1-δ) methods. Moreover, in this case the alternate 2-δ method also has the
	same number of iterations as the alternate secant method (see Table 1) even if

Table 4 .

 4 Again, no numerical order of convergence can be clearly drawn from this test case. crossed δ 2 alternate δ 2 crossed 2-δ alternate 2-δ crossed 2-δ bis

	ε=10 -6	431	287	326	129	Not Cvg
	ε=10 -8	940	678	579	204	Not Cvg
	ε=10 -12	Not Cvg	1691	1613	476	Not Cvg

Table 4 :

 4 Total number of iterations involved for the Polycrystals test case -Three iteration residual methods

In practice, other minimization strategies are often used such as the Gram-Schmidt orthogonalization algorithm or the QR decomposition. This latter method seems in this case to be the most efficient[START_REF] Degroote | Vierendeels, Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction[END_REF][START_REF] Walker | Anderson acceleration for fixed-point iterations[END_REF].

Acknowledgements

This work was accomplished within the framework of the PLEIADES project, financially supported by the CEA (Commissariat à l'Énergie Atomique et aux Énergies Alternatives), EDF (Électricité de France) and AREVA.

The test cases used in Section 5 were implemented by Jean-Michel Proix.

Authors are grateful to P. Verpeaux and the whole Cast3M development team for fruitful discussions on their use of acceleration methods in the context of structural mechanics.

is lower than a prescribed accuracy ε. This stopping criterion guarantees that the obtained solution is the fixed point equation solution.

MTest can be used to simulate simple mechanical tests. In this case, ε is usually chosen equal to 10 -6 or 10 -8 to obtain satisfactory results. This corresponds to a relative precision lower than 10 -3 or 10 -5 as the strains typically have a magnitude ranging between 10 -3 and 10 -1 .

A practical way of checking that the behaviour integration gives repeatable results whatever compilers, compiler optimisation flags and/or operating systems is to use a ε value smaller than the criterion values used within the behaviour integration. Typically, a stringent value of 10 -12 , corresponding to a relative precision of about 10 -9 , is used for ε (precision by default) while the behaviour integration is performed with a criterion of 10 -8 .

Test cases

In order to illustrate the performances of the various algorithms described in this paper, we have chosen two test cases, called Chaboche and Polycrystals, respectively and which shall now be briefly described.

Chaboche test case

This test case describes an isotropic standard plastic behaviour with two kinematic hardening rules derived from the works of J.L. Chaboche [START_REF] Chaboche | Constitutive equations for cyclic plasticity and cyclic viscoplasticity[END_REF][START_REF] Chaboche | Cyclic viscoplastic constitutive equations[END_REF]. Such a behaviour is commonly used to describe metal plasticity in engineering.

A whole description of the behaviour can be found in the Code-Aster documentation [START_REF] Edf | Elasto-visco-plastic Chaboche constitutive law, Code-Aster reference guide R5.03.04[END_REF].

The behaviour is time-independent, so the time scale is arbitrary. The test case imposes two components of the stresses: σ xx and σ xy which monotonically and linearly increase from 0 at the beginning of the test to 143.5 MPa at the end of the test. The mechanical equilibrium imposes that the other components of the stresses are null. This loading is described in 13 time steps. The results of the simulation are reported in Figure 1.

The residual methods with M = 1 and Z 1 n = δ 2 X n-i-1 as well as the residual methods with M = 2 and Z i n = δX n-i , i = 1, 2 are based on the last three iterations and will also be compared:

• the alternate δ 2 method

• the crossed 2-δ method

• the alternate 2-δ method (equivalent to the Anderson extrapolation method with M = 2 [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF])

It can be easily shown that the alternate 2-δ method remains the same if we consider G(X n ) -G(X n-2 ) instead of G(X n-1 ) -G(X n-2 ) in the expression of X n+1 . The new coefficients λ1 n and λ2 n are hence linear combinations of λ 1 n and λ 2 n :


Abstract :

The Pellet-Cladding mechanical Interaction (PCI) occurs in Pressurized Water Reactors and might pose a risk to the integrity of the cladding containing the fuel. In this paper, the Local Defect Correction multigrid method (LDC) is performed in order to improve the PCI simulation accuracy. From an initial coarse mesh, this method consists in recursively adding local sub-grids in zones where higher accuracy is required. In order to automatically detect the zones to be refined, the LDC method is coupled with the Zienkiewicz-Zhu a posteriori error estimator. This estimator is based on the stress discontinuities between elements obtained by the finite element method. It provides hence a tool to detect the stress singularity zones. The strategy proposed enables us to automatically refine the mesh around the stress concentration zones and to stop the refinement process when the solution is accurate enough or when the refinement does not improve the solution accuracy anymore. The numerical results on 2D elasticity test cases show the efficiency of this strategy but also the limit of the Zienkiewicz-Zhu a posteriori error estimator.

Mots clefs : Automatic mesh refinement, Local Defect Correction method, Zienkiewicz and Zhu a posteriori error estimator, Stress concentration localization 1 Introduction

In Pressurized Water Reactors, the Pellet-Cladding mechanical Interaction (PCI) happens in early stage of irradiation. Precise numerical simulations of this localized phenomenon require cells of 1 µm for a structure of 1 cm. Actual personal computers are not able to reach this accuracy in reasonable computational time and memory space. In order to overcome these difficulties, it is proposed here to use the Local Defect Correction (LDC) multigrid method introduced by Hackbush [START_REF] Hackbusch | Local Defect Correction Method and Domain Decomposition Techniques[END_REF]. This method consists in recursively adding local sub-grids in zones where higher accuracy is required. In order to automatically detect the zones to be refined, a coupling between the LDC multigrid method and the Zienkiewicz-Zhu (ZZ) a posteriori error estimator [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF] is proposed. A strategy for coupling the LDC multigrid method with the ZZ estimator was already introduced in [START_REF] Barbié | Strategies around the local defect correction multi-level refinement method for three-dimensional linear elastic problems[END_REF][START_REF] Barbié | An automatic multilevel refinement technique based on nested local meshes for nonlinear mechanics[END_REF]. However, strategy introduces an arbitrary coefficient which may depend on the problem under consideration and on the user knowledge . Another strategy, proposed in [START_REF] Liu | Local defect correction method coupled with the zienkiewiczzhu a posteriori error estimator in elastostatics solid mechanics[END_REF] and recalled in this paper, has been developed in order to avoid this dependence. The numerical results presented in [START_REF] Liu | Local defect correction method coupled with the zienkiewiczzhu a posteriori error estimator in elastostatics solid mechanics[END_REF] on 2D axisymmetric test cases show that this method works well. The numerical results for the 2D plane strain test case presented in this paper confirm the efficiency of this automatic strategy but also show the limit of this method due to the ZZ error estimator.

This paper is organized as follows. The section 2 focuses on the coupling of the LDC method with the ZZ a posteriori error estimator. The section 3 introduces briefly the PCI phenomenon and the associated simplified 2D plane strain test case. Finally, numerical results that demonstrate the automatic detection and refinement of stress concentration zones are presented in section 4.

2 Automatic mesh refinement strategy

The LDC method

The Local Defect Correction method (LDC) [START_REF] Hackbusch | Local Defect Correction Method and Domain Decomposition Techniques[END_REF] method is an Adaptive Mesh Refinement method based on a local multigrid process [START_REF] Angot | Une méthode adaptative de raffinement local : la Correction de Flux à l'Interface[END_REF][START_REF] Brandt | Multi-level adaptive solutions to boundary-value problems[END_REF]. This methods starts with a global coarse grid defined on the whole domain. Then local sub-grids with finer and finer mesh steps are recursively added in zones of interest. Prolongation and restriction operators are used to link the different levels of grids. The problems arising on all the grids are solved sequentially until the solution of the coarsest grid has converged. This iterative process can be described by a ∧-cycle, see Figure 1. In the following, the two-grids LDC method is briefly presented. The reader can refer to [START_REF] Hackbusch | Local Defect Correction Method and Domain Decomposition Techniques[END_REF][START_REF] Barbié | Strategies around the local defect correction multi-level refinement method for three-dimensional linear elastic problems[END_REF] for more details.

A problem (P) well-defined on an open-bounded domain Ω with boundary ∂Ω is considered. This problem is discretized on a coarse grid G H of boundary Γ H with grid size H. The discretized problem is denoted (P H ). We define a local subdomain Ω 1 ⊂ Ω covering the zone of interest of Ω. Ω 1 is discretized with a local grid G h of boundary Γ h with grid size h < H (see Figure 2). The discretization of (P) on G h is denoted (P h ).

The LDC algorithm then writes :

-Initialization : Solve the discretized problem (P H ) on G H .

-Until the solution of (P H ) converges between two iteration cycles do : -Prolongation step : Define the boundary conditions of the local problem (P h ). On the boundary part The multigrid algorithm is obtained by recursively applying the two-grids algorithm [START_REF] Hackbusch | Local Defect Correction Method and Domain Decomposition Techniques[END_REF][START_REF] Barbié | Strategies around the local defect correction multi-level refinement method for three-dimensional linear elastic problems[END_REF]. This method has the main advantage to consider the solvers as "black-box". It should also be noted that it is a generic method : solver, refinement ratio, mesh type, model, etc, could be different (or not) between the different levels of refinement. However, with the definition of the correction zone, this method requires a sufficiently large refinement zone to be efficient [START_REF] Hackbusch | Local Defect Correction Method and Domain Decomposition Techniques[END_REF].

Γ h ∩ ∂Ω,

Zienkiewicz and Zhu a posteriori error estimators

A posteriori error estimators [START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF][START_REF] Ladevèze | Error estimate procedure in the finite-element method and applications[END_REF][START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF] provide a technique to automatically detect the elements with higher error from a preliminary calculation.

Based on the fact that the stress field of the Lagrange finite element method is discontinuous between elements, recovery-based error estimators proposed by Zienkiewicz and Zhu [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimation. Part I : The recovery technique[END_REF] intend to construct a smoothed stress field. These estimators propose to recover the nodal values of stresses by using a minimization method and then construct a continuous stress field using displacement shape functions. The error estimator is evaluated by the difference between the finite element stress solution and the smoothed stress field. Because of its efficiency, the good results it provides and the easiness of its implementation, these estimators are widely used in industrial codes such as Cast3M or Code_Aster.

Coupling the LDC with the ZZ estimators

In order to automatically detect the zones of interest, we couple the ZZ estimator within the LDC method. The coupling between the ZZ estimator and the LDC method can be easily realized as both approaches are sequentially applied : at the first prolongation step, after each calculation, the ZZ estimator is applied on each level to detect the zones of interest. These sub-grids are then used in the sequel of the LDC process.

In order to automatically detect the zones of interest with the error indicator estimated by ZZ estimator, Barbie et al [START_REF] Barbié | Strategies around the local defect correction multi-level refinement method for three-dimensional linear elastic problems[END_REF][START_REF] Barbié | An automatic multilevel refinement technique based on nested local meshes for nonlinear mechanics[END_REF] have proposed to refine the elements K for which the indicator is superior to α of the maximum error :

e K > α(max L⊂G l e L ) (1)
with 0 ≤ α ≤ 1, α a unknown coefficient
This strategy is easy to implement but requires the knowledge of the factor α which may depend on the considered problem. Moreover, an arbitrary stopping criteria has to be defined because this indicator never stops.

The following automatic refinement procedure is proposed (see also [START_REF] Liu | Local defect correction method coupled with the zienkiewiczzhu a posteriori error estimator in elastostatics solid mechanics[END_REF]) to overcome these disadvantages.

The error between the finite element solution and the smoothed stress fields used by the Zienkiewicz and Zhu a posteriori error estimators calculated in energy norm writes :

e K E = K (σ * -σ h ) : (ε * -h ) dK K σ * : ε * dK 1/2 (2)
with σ * smoothed stress solution, σ h the finite element stress solution, C is the stiffness tensor and

ε * = [C] -1 σ * , ε h = [C] -1 σ h .
The automatic procedure of refinement proposed in [START_REF] Liu | Local defect correction method coupled with the zienkiewiczzhu a posteriori error estimator in elastostatics solid mechanics[END_REF] works directly on the elements for which the stress error is superior to a relative threshold β :

e K E > β (3)
When the refinement stops, the estimated global error e Ω E is guaranteed to be inferior to the relative threshold β.

When the problem presents a singularity, the error between the finite elements solution and the smoothed solution could be always superior to the relative threshold. So an additional stopping criteria for the automatic generation of the sub-grids has to be defined : when the measure of the detected zone is inferior to the measure of one coarsest element, the generation of sub-grids is stopped. In fact, following the LDC algorithm (see section 2.1), the local zone detected in this case is too small to correct the coarsest solution.

3 The PCI phenomenon and the associated simplified test case

The PCI phenomenon

The PCI phenomenon occurs during the irradiation in Pressurized Water Reactors. The fuel pellets crack at the early stages of irradiation. Meanwhile, the pressure of the water in the primary circuit leads to the cladding's creep and the pellets swell due to the accumulation of gaseous fission products within the material. The initial gap between the pellets and the cladding decreases and discontinuous contacts appear. At the same time, the fragmented pellets deform in a hourglass shape due to high temperature gradient. All these phenomena induce localized high stress concentrations on the cladding in front of the inter-pellet plane that might lead to fuel rod failure under certain circumstances.

The simplified 2D(r,θ) test case

A test case with the plane strain hypothesis, called the 2D(r,θ) test case, is introduced. This test case is designed to simulate the effect of the pellet cracking on the cladding. Only the cladding is considered here and the contact with the pellet is represented by a pressure discontinuity (due to the cracks) on the internal radius of the cladding. We focus on the elastic response of the cladding with the Young's modulus E = 100 GP a and the Poisson's ratio ν = 0.3. As the pellet is assumed to crack in a regular way (see figure 5a), only 1/16 of the cladding is modelled. The definition of the problem is represented on Figure 5b. A reference solution obtained with a very fine mesh (quasi-uniform mesh of 2µm mesh step) is used to validate the proposed strategy.

Numerical results

We apply the strategy of coupling LDC method with ZZ estimator introduced in the section 2.3 to the 2D(r,θ) test case. Figure 6 presents the generated sub-grids for an initial mesh size equal to 54.5 µm and a threshold β of 1%. As expected, the local sub-grids are located around the pressure discontinuity zone where the stresses present a singularity.

This example shows that the proposed coupling strategy works very well in automatically detecting the zones to be refined. Then the LDC algorithm is applied on these sub-grids. The iterative process is considered as convergent when the L 2 norm relative error between two successive coarse solutions u 0 is smaller than 1.10 -5 :

u k 0 -u k-1 0 L 2 u k 0 L 2 ≤ 1.10 -5 (4)
The LDC method converges very fast. Typically, less than four iterations are required for convergence.

The coupling strategy has been tested on different coarse meshes and different thresholds. The results are reported in Table 1. The results of each calculation are presented in three lines : the first line indicates the number of sub-grid automatically generated, the second line contains the numbers of nodes of each subgrid while the third line gives the relative error in energy norm between the convergent LDC composite solution and the reference solution. The first conclusion to be drawn from Table 1 is that the generation of sub-grids automatically stops. Table 1 also shows the efficiency of the LDC method. The errors are quickly reduced by adding local grids with a limited number of degrees of freedom. So it is more interesting to begin with a coarse enough grid then add local sub-grids rather than begin with a too fine mesh. Moreover, at the end of the calculation, the errors obtained are generally lower than the thresholds. Concerning the ZZ error estimator, it seems that for coarse meshes (see row 1 of Table 1) the estimator underestimates the finite element error while for fine meshes this estimator overestimates the finite element error (see row 3

of Table 1 where more sub-grids are generated with the decrease of the threshold whereas it was not necessary).

However, only the results for the initial coarsest mesh (with mesh size equal to 109 µm) are not as good as expected with the three largest thresholds (5%, 2%, 1%). To find out the reason of this behaviour, the relative error estimated by the ZZ estimator and the relative stress error between the finite element solution and the reference solution are compared on Figure 7. For this coarsest mesh, the major error between the finite element solution and the reference solution is located around the bottom right corner of the mesh. This is due to a bad approximation of the boundary condition. This kind of error could not be detected by the ZZ estimator which focuses only on the smoothness of the solution. In these zones, the ZZ estimator underestimates the finite element error. This figure can also explain why the result for the same initial mesh size with a threshold β of 0.5% is much better. Indeed, in this case, the strategy leads to refine all the domain and hence the boundary condition is better approximated.

Conclusion

In order to accurately simulate the stress concentration zone of the PCI, a strategy to couple the Local Defect Correction method with the Zienkiewicz and Zhu a posteriori error estimator in an user-independent way has been proposed. This strategy works well for the test cases whose errors are mainly caused by a stress singularity. The LDC method appears to be a powerful adaptive mesh refinement tool. The ZZ error estimator well detects stress singularity zones but is unable to detect the zones where the smooth solution is not accurate enough. For these situations, it could be interesting to combine another error estimator with the LDC method in order to more accurately detect all the zone to be refined.

F 5 -

 5 2D (r,θ) test case

F 7 -

 7 Error estimated by the ZZ estimator (at left) and error between the finite element solution and the reference solution (at right) (mesh size : 109 µm)

Automatic multigrid mesh refinement strategy for precise localized stress concentration simulation H. LIU a,b , I. RAMIÈRE a , F. LEBON b a

  . DEN/DEC/SESC/LSC, CEA Cadarache, F-13108 Saint-Paul Lez Durance, France b. Laboratoire de Mécanique et d'Acoustique, CNRS, UPR 7051, Aix-Marseille Univ, Centrale Marseille, 31, Chemin Joseph Aiguier, F-13402 Marseille Cedex 20, France

	Résumé :

L'interaction mécanique pastille-gaine est un phénomène se produisant dans les réacteurs à eau pressurisée et pouvant mettre à défaut l'intégrité de la gaine entourant le combustible. Dans cet article, la méthode multigrille Local Defect Correction (LDC) est utilisée afin d'améliorer la précision de la simulation de ce phénomène. A partir d'un maillage initial, cette méthode consiste à ajouter récursivement des sous-grilles locales dans des zones où une solution précise est recherchée. Afin de détecter automatiquement les zones à raffiner, la méthode LDC est couplée avec l'estimateur d'erreur a posteriori de type Zienkiewicz et Zhu. Cet estimateur se base sur le fait que les contraintes obtenues par la méthode des éléments finis sont discontinues entre éléments. Il permet ainsi de détecter les zones de singularité de contraintes. La stratégie proposée permet de raffiner automatiquement le maillage dans les zones de concentration de contraintes et également d'arrêter le raffinement quand la solution est suffisamment précise ou lorsque le raffinement n'apporte plus d'amélioration de la solution. Les résultats numériques obtenus sur des cas tests 2D élastiques avec discontinuité de chargement montrent l'efficacité de cette stratégie mais également les limites de l'estimateur d'erreur a posteriori de type

Zienkiewicz et Zhu.

  Solve the local discretized problem (P h ). -Restriction step : Calculate the defect in order to correct the coarse problem (P H ). The fine problem (P h ) solution is restricted on the set A H of coarse nodes of G H strictly included in G h . The defect of the problem (P H ) is then evaluated on the interior nodes of A H ∪ Γ H in the sense of the discretization scheme, see Figure 4. -(P H ) is then solved by adding the defect to the initial source term.

					Fine grid G l *
						Fine grid G	h
					Coarse grid G 0
	Initialisation		Prolongation step (boundary conditions)
	Smoothing or exact solving		Restriction step (correction)	Coarse grid G	H	Zone of interest
	Converged solution				F	2 -Example of hie-
						rarchical meshes used in LDC
	F	1 -Representation of the iterative process in ∧-cycle	method
	Original problem boundary conditions -Porjection of next coarser solution	1
	F	3 -Example of boun-	F	4 -Example of restriction zone A H (on the left) and
	dary conditions on G h for hie-	correction zone (on the right) in the case of a 5-point stencil
	rarchical meshes		operator for hierarchical meshes

original problem (P) boundary conditions are imposed. On Γ h \∂Ω, Dirichlet boundary conditions are obtained by projecting the coarse problem (P H ) solution, see Figure

3

.