Christophe Bovet

Augustin Parret-Fréaud

Pierre Gosselet

Two-level adaptation for Adaptive Multipreconditioned FETI

Keywords: Domain decomposition, FETI, Krylov subspace methods, multipreconditioning, adaptive multipreconditioning, AMPFETI, SFETI See [6]

This article introduces two strategies to reduce the memory cost of the Adaptive Multipreconditioned FETI method (AMPFETI) while preserving its capability to solve ill conditioned systems efficiently. Their common principle is to gather search directions into aggregates which are frequently adapted in order to achieve the best compromise between the decrease of the solver error and the computational resources employed. The methods are assessed on two weak scalability studies on highly heterogeneous problems up to 10368 cores and half a billion of unknowns, and on two illconditioned industrial applications, related to the numerical homogenization of solid propellant and to the simulation of a multiperforated aircraft combustion chamber.

Introduction

The continuous increase of high performance computing resources naturally appeals the engineers to take into account more complex physical phenomena and to simulate fine scale models. In the area of solid mechanics, for instance, this can be exemplified by the inclusion of small geometric details (microperforations, convoluted blade cooling channels. . . ), or complex nonlinear heterogeneous materials (woven composite, crystal plasticity. . . ). After the discretization of the involved partial differential equations, this results in large sparse linear(ized) systems which range from millions to billions of unknowns, and which are very often ill-conditioned.

Non overlapping domain decomposition methods, such as the Balancing Domain Decomposition (BDD) [START_REF] Mandel | Balancing domain decomposition[END_REF] or the Finite Element Tearing and Interconnecting (FETI) [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] and their constrained variants FETI-DP [START_REF] Farhat | FETI-DP: a Dual-Primal Unified FETI Method -Part I: a Faster Alternative to the Two-Level FETI Method[END_REF] and BDDC [START_REF] Dohrmann | A preconditionner for substructuring based on constrained energy minimization[END_REF], provide a favorable framework to solve these linear systems in parallel. Recently developed domain decomposition methods, such as the ML-FETI-DP [START_REF] Toivanen | A multilevel FETI-DP method and its performance for problems with billions of degrees of freedom[END_REF] or the HTFETI [START_REF] Merta | Intel xeon phi acceleration of hybrid total feti solver[END_REF], are able to tackle really large-size problems while exhibiting a really good scalability. Indeed, these methods are well adapted to supercomputers since they rely on both iterative and direct solvers. Usually each subdomain is assigned to one computer node (or one socket) where a multithreaded direct solver is used. An interface problem is solved between subdomains with a Krylov solver (MPI parallelism) to obtain the global solution.

Thanks to this combination of direct and iterative solvers, domain decomposition methods are usually more robust than traditional iterative solvers as they ensure a logarithmic bound in H/h on the condition number, where h is the characteristic size of the finite element and H that of a subdomain [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF]. Despite that, real engineering applications frequently exhibit pathological components that hinder the convergence of the underlying Krylov solver, such as jagged interfaces, bad aspect ratios, strong heterogeneities misplaced with respect to the interface, incompressibility etc.

Two kinds of strategy have been proposed to increase the robustness of domain decomposition solvers. The first strategy is to detect a priori the local contributions which are known to penalize the convergence, and to remove them by augmenting the Krylov solver. The problematic local information can be identified through local generalized eigenvalue problems [START_REF] Spillane | Automatic spectral coarse spaces for robust FETI and BDD algorithms[END_REF][START_REF] Spillane | Robust domain decomposition methods for symmetric positive definite problems[END_REF][START_REF] Mandel | Adaptive selection of face coarse degrees of freedom in the BDDC and the FETI-DP iterative substructuring methods[END_REF] or reused from a previous solution [START_REF]Conjugate Gradient Method with Preconditioning by Projector[END_REF][START_REF] Farhat | Extending substructure based iterative solvers to multiple load and repeated analyses[END_REF][START_REF] Risler | Iterative accelerating algorithms with Krylov subspaces for the solution to large-scale non-linear problems[END_REF][START_REF] Gosselet | Total and selective reuse of Krylov subspaces for the solution to a sequence of nonlinear structural problems[END_REF][START_REF] Michael | Recycling of Solution Spaces in Multi-Preconditioned FETI Methods Applied to Structural Dynamics[END_REF]. The second strategy is to use a block Krylov solver [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF][START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF] or a multipreconditioned Krylov solver [START_REF] Bridson | A multipreconditioned conjugate gradient algorithm[END_REF][START_REF] Spillane | An Adaptive Multipreconditioned Conjugate Gradient Algorithm[END_REF][START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF]. Multipreconditioning is well adapted to the FETI and BDD methods because of the additive nature of their classical preconditioners. Moreover, this technique has been successfully applied to solve non symmetric systems [START_REF] Greif | GMRES with multiple preconditioners[END_REF][START_REF] Bovet | Multipreconditioning for nonsymmetric problems: The case of orthomin and biCG[END_REF] for which finding the suitable augmentation is still an open question.

The fast growing search space produced by a multipreconditioned solver is the key ingredient to ensure a fast convergence. In the original Simultaneous FETI method (SFETI, which corresponds to multipreconditioned conjugate gradient applied to FETI) [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF], the algorithm generates at each iteration as many search directions as there are subdomains in the decomposition. However, these extra search directions need to be stored in memory which can become problematic in case of a large number of subdomains. To reduce the memory requirement of SFETI, a sparse storage format has recently been proposed in [START_REF] Molina | New implementations for the Simultaneous-FETI method[END_REF]. An alternative trail is followed in Adaptive Multipreconditioned FETI (AMPFETI) [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF][START_REF] Bovet | Multipreconditioning for nonsymmetric problems: The case of orthomin and biCG[END_REF][START_REF] Bovet | Méthode de décomposition de domaine multipréconditionnée et adaptative pour les problèmes mal conditionnés[END_REF]: among the search directions generated by the subdomains, only the ones that are identified as useful for the convergence are kept individually, the others are reduced by summation as done in classical FETI [START_REF] Spillane | An Adaptive Multipreconditioned Conjugate Gradient Algorithm[END_REF]. Thanks to this test, the size of the generated search space is significantly reduced compared to SFETI. In [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF], the AMPFETI's capability to solve ill-conditioned systems has been assessed up to 500 subdomains and 100 millions of unknowns. Note that for now, multipreconditioned solvers have been tested with FETI and BDD domain decomposition methods, but they could also be employed with the FETI-DP and BDDC methods.

The aim of this paper is to present two improvements of AMPFETI allowing addressing larger illconditioned problems. First, a better control of the memory requirement is achieved by aggregating search direction candidates before multipreconditioning. The paper provides guidelines for the construction of the aggregates. Second, in order to accelerate the convergence even more, aggregates are updated at each iteration in order to generate better search directions.

The article is organized as follows: Section 2 gives a short presentation of Adaptive Multipreconditioned FETI; then Section 3 details the methodology proposed to gather the search directions into adaptive aggregates. Scalability results on academic benchmarks are provided in Section 4 and Section 5 presents two realistic industrial simulations. Section 6 concludes the paper.

Adaptive multipreconditioned FETI: a short reminder

Adaptive multipreconditioned FETI is the combination of the FETI domain decomposition method [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | The dual Schur complement method with well-posed local Neumann problems[END_REF] and the adaptive multipreconditioned conjugate gradient solver [START_REF] Spillane | An Adaptive Multipreconditioned Conjugate Gradient Algorithm[END_REF][START_REF] Bridson | A multipreconditioned conjugate gradient algorithm[END_REF].

The FETI linear system

Let us consider the linear system of equation Ku = f arising from the finite element discretization of a linear mechanical problem defined on domain Ω. Let u denote the vector of generalized displacement and f the vector of generalized forces. In the present study, we assume that the stiffness matrix K is symmetric positive definite. Let (Ω s ) 1 s N d be a non overlapping partition of Ω such that: ∀s = p, Ω s Ω p = ∅ and Ω =

N d s=1

Ωs where Ωs is the closure of Ω s . The interface between two subdomains s and p is denoted by Υ sp = Ωs Ωp and the union of all interfaces of the subdomain s is denoted by Υ s . The set gathering all interfaces is denoted by Υ.

In the FETI method, a Lagrange multiplier field λ that connects subdomains is introduced and the global linear system Ku = f is substituted by an equivalent substructured formulation

K s u s =f s + T s B s λ s ∀ 1 s N d (1) N d s=1 B s T s u s =0 (2)
where T s : Ω s → Υ s are trace operators and B s are signed Boolean assembly operators on the global interface. From a mechanical point of view, equations (1) express the equilibrium of all subdomains and (2) corresponds to the continuity of the displacement across interfaces. Then, the local equilibrium is expressed in terms of interface unknown. Schur complements (S s ) naturally arise and, after some algebraic manipulations, the classical FETI system is formed:

F G G 0 λ β = d e (3)
where

e = -f 1 R 1 | . . . |f N d R N d ; d = - N d s=1 B s T s K s † f s F = N d s=1 B s S s † B s ; S s = K s bb -K s bi K s -1 ii K s ib G = B 1 T 1 R 1 | . . . |B N d T N d R N d
The Schur complement S s is computed by condensing out the internal degrees of freedom (index i) on the boundary (index b). S s † denotes the generalized inverse of S s in case the kernel R s of K s is not restricted to zero. The second row, G λ = e, is associated with the constraint that each subdomain must remain self-equilibrated. The vector β contains the rigid body mode coefficients of each subdomains. This vector is computed after λ has been found (see [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]). Please note that the FETI operator F is never assembled since only the result of a multiplication by F is needed. Instead of solving a saddle point system, an initialization-projection strategy is applied and λ is sought as:

λ = λ 0 + Π λ ; G λ 0 = e ; G Π = 0
The projector Π is of the form Π = I -AG(G AG) -1 G , where A is a symmetric positive definite matrix most often chosen among the classical preconditioners S and the identity matrix.

Substituting this form into (3), and pre-multiplying by Π leads to the final linear system:

Π F Π λ = Π (d -F λ 0 ) (4)
In the classical FETI method, the system (4), is solved with a preconditioned conjugate gradient. The usual preconditioner is a scaled sum of (approximation of) Schur complements:

S = N d s=1
Bs Šs Bs [START_REF] Bovet | Méthode de décomposition de domaine multipréconditionnée et adaptative pour les problèmes mal conditionnés[END_REF] where ( Bs ) are scaled assembly operators such that:

N d s=1 B s Bs B j = B j , ∀ 1 j N d
Local approximations Šs of S s can be chosen among S s (Dirichlet preconditioner), K s bb (Lumped preconditioner) and diag(K s bb ) (Superlumped preconditioner). In the following, we will single out the subdomains' contribution to the operator and to the preconditioner, and we will use the notation:

F s = B s S s † B s ; Ss = Bs Šs Bs

Multipreconditioned FETI

There are many cases in which the preconditioned conjugate gradient used in the FETI algorithm converges very slowly or does not converge at all, even with state-of-the-art choices of the preconditioner. A typical pathological case is strong material heterogeneity misplaced with respect to the interface between subdomains [START_REF] Daniel | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF]. Such a situation cannot be avoided for the simulation of solid propellant (see Section 5.1).

Recently, multipreconditioned strategies [START_REF] Rixen | Substructuring and Dual Methods in Structural Analysis[END_REF][START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF][START_REF] Bovet | Méthode de décomposition de domaine multipréconditionnée et adaptative pour les problèmes mal conditionnés[END_REF][START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF] were proposed to solve these hard problems. In these methods, the contribution of each subdomain to (5) is considered as a separate preconditioner which, at each iteration, supplies one direction inside a search space of dimension N d . All these strategies rely on multipreconditioned Krylov solvers [START_REF] Greif | GMRES with multiple preconditioners[END_REF][START_REF] Bovet | Multipreconditioning for nonsymmetric problems: The case of orthomin and biCG[END_REF] and more specifically multipreconditioned conjugate gradient (MPCG) [START_REF] Bridson | A multipreconditioned conjugate gradient algorithm[END_REF].

If r i+1 denotes the residual at iteration i+1, the classical preconditioned conjugate gradient generates a search space of dimension 1 materialized by the preconditioned residual z i+1 , multipreconditioned conjugate gradient generates a larger search space, spanned by N d vectors gathered in the matrix Z i+1 :

z i+1 = N d s=1 Ss r i+1 ; Z i+1 = S1 r i+1 , . . . , SN d r i+1 (6)
The optimality conditions of conjugate gradient enables to compute the best combinations of the columns of Z i+1 at each iteration, whereas classical preconditioned conjugate gradient corresponds to the fixed choice z i+1 = Z i+1 1 (with 1 = 1, . . . , 1 ).

Algorithm 1 presents multipreconditioned FETI combined with the adaptive τ -test and an aggregation procedure described in Section 3. To recover the original multipreconditioned FETI, forget the τ -test of line 12, consider In [START_REF] Michael | Recycling of Solution Spaces in Multi-Preconditioned FETI Methods Applied to Structural Dynamics[END_REF] it was shown that multipreconditioning permits to explore the subspace spanned by the "bad eigenvectors" detected by the Generalized Eigenvalues in the Overlaps (GenEO) algorithm, and thus it leads to a more stable convergence on ill-conditioned problems. But multipreconditioning comes with a high price. Indeed, the original multipreconditioned FETI generates at each iteration, two matrices W i and Q i (see Algorithm 1). Since, a full reorthogonalization is mandatory, all these matrices need to be stored.

C k = {Ω k } ∀ 1 k N a = N d ,
To set the ideas, let us consider a regular cubic mesh with a cubic domain decomposition, and let n denote the number of nodes per subdomain edge (see Figure 1 and Figure 3). The memory footprint of all the elements of the sequences

(W i ) i and (Q i ) i is O(n 2 N d n iter ) and for (∆ i ) i it is O(N 2
d n iter ). This may become prohibitive when the number of subdomains N d scales with the same rate as n does.

As mentioned in the introduction, one possibility to reduce the memory requirement of MPFETI is to express all operations in terms of (Z i ) i instead of (W i , Q i ) i . Indeed, contrarily to (W i , Q i ) i , (Z i ) i is sparse (the contribution of one subdomain only affects its neighbors). (Z i ) i is thus much less expensive to store. This strategy, described in [START_REF] Molina | New implementations for the Simultaneous-FETI method[END_REF], comes with many algebraic manipulations.

Algorithm 1: AMPFETI with local τ -test and aggregation

Input: Convergence threshold > 0, adaptation threshold τ > 0, initial partition

PN a 1 Π = I -AG(G AG) -1 G 2 λ0 = AG(G AG) -1 e 3 r0 = Π (d -F λ0) 4 Z0 = Ω s ∈C 1 Ss r0| . . . | Ω s ∈C Na Ss r0 5 W0 = ΠZ0 6 λ0 = 0, i = 0 7 while ri > do 8 Qi = F Wi 9 ∆i = Q i Wi ; γ i = Z i ri ; αi = ∆ † i γ i 10 λi+1 = λi + Wiαi 11 ri+1 = ri -Π Qiαi 12 t C i = (Wiαi) Ω s ∈C F s (Wiαi) r i+1 Ω s ∈C Ss ri+1 13 Zi+1 = concatenate({ Ω s ∈C Ss ri+1 ; t C i < τ }, t C i >τ Ω s ∈C Ss ri+1) 14 Wi+1 = ΠZi+1 for 0 ≤ j ≤ i do Φi,j = Q j Wi+1
Wi+1 ← Wi+1 -Wj∆ † j Φi,j 15 return λ0 + Πλi

Adaptive multipreconditioned FETI

A major improvement of MPCG is the addition of an adaptive process [START_REF] Spillane | An Adaptive Multipreconditioned Conjugate Gradient Algorithm[END_REF] which allows singling out useful search directions. Before being inserted in the search subspace, a criterion enables to forecast the effectiveness of the search directions provided by the subdomains with regard to the convergence. Depending on a user parameter τ which can be linked with a target relative reduction of the F -norm of the error, the adaptive criterion chooses between 1 and N d search directions. In [START_REF] Spillane | An Adaptive Multipreconditioned Conjugate Gradient Algorithm[END_REF], two variants of the adaptive process were proposed: The adaptive algorithms [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF] start from the full subspace generated by multipreconditioning [START_REF] Bovet | Two-level adaptation for adaptive multipreconditioned FETI[END_REF], but now the columns are only considered as potential search directions (hence the superscript p):

Z p i+1 = S1 r i+1 , . . . , SN d r i+1 (7)
With the Global τ -test, the effect of multipreconditioning can be forecast by the quantity (with Algorithm 1 notations):

t G i = γ i α i r i+1 z i+1 (8)
which is almost costless to compute. If t G i is less than a certain criterion τ then it is worth multipreconditioning using Z p i+1 as search space, else the classical FETI approach with one dimensional search space z i+1 already performs satisfyingly.

With the Local τ -test, the contribution to the decrease of the error of each vector resulting from the local preconditioning Ss r i+1 can be estimated individually. This estimate is given by the scalar t s i :

t s i = (W i α i ) F s (W i α i ) r i+1 Ss r i+1 ( 9
)
where all quantities defined at iteration i are known. Note that computing the local test is significantly more expensive than the global one. If t s i is less than a certain criterion τ , then the vector is worth being isolated, else it can be harmlessly combined with other little-contributing vectors. Let I i be the set of subdomains (Ω s ) such that t s i < τ , the actual search subspace then can be written as:

Z i+1 =   . . . , Ss r i+1 , . . . s∈Ii , s / ∈Ii Ss r i+1   (10)
We refer to [START_REF] Spillane | An Adaptive Multipreconditioned Conjugate Gradient Algorithm[END_REF] and [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF] to get more details about the adaptive tests. Algorithm 1 presents Local AMPFETI combined with an aggregation procedure which will be described in Section 3. To recover the original Local AMPFETI, just consider

C k = {Ω k } ∀ 1 k N a = N d .
It is worth noting that unfortunately, AMPFETI cannot be combined with the sparse storage strategy of [START_REF] Molina | New implementations for the Simultaneous-FETI method[END_REF]. Indeed, for a given iteration j, all the columns (before selection) of the previous preconditioned residual (. . . , Ss r i , . . .) s;i<j need to be stored in order to rebuild the selected search directions, even if they are scarce.

Two-level adaptive strategy 3.1 Gathering search direction into aggregates

Even with a very stringent test, the memory requirement of the original AMPFETI may become high when the system to be solved is ill-conditioned and the number of subdomains is large.

As a complement to the selection threshold τ , we propose to choose a priori the granularity of the multipreconditioning by gathering search directions into aggregates. Indeed, the MPCG framework is not limited to starting from individual contributions from subdomains as in [START_REF] Bridson | A multipreconditioned conjugate gradient algorithm[END_REF].

Let P Na = (C k ) 1≤k≤Na be a partition of the subdomains' connectivity graph into N a N d aggregates. We can then propose N a potential search directions under the form:

Z p i+1 =   Ω s ∈C1
Ss r i+1 , . . . ,

Ω s ∈C Na Ss r i+1   (11)
Then, the adaptive τ -test can be applied on each aggregate's contribution, as described in Algorithm 1.

In the end, at most N a vectors are retained and used by the solver, then the memory footprint for the storage of (W i ) i and An interesting point is that the one-aggregate case (N a = 1) corresponds to the genuine FETI method equipped with a classical conjugate gradient.

(Q i ) i can not exceed O(n 2 N a n iter ) and O(N 2 a n iter ) for (∆ i ) i .
In practice, the partition (C k ) 1≤k≤Na is provided by a graph partitioning library such as Metis [START_REF] Karypis | A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs[END_REF]. More precisely, let G = Ω 1 , . . . , Ω N d , Υ be the subdomain's connectivity graph. Subdomains (Ω s ) are the graph vertexes. There is an edge between two subdomains (Ω s , Ω p ) if they share one interface Υ sp .

The choice of optimal weights is hard to now a priori. In this Section, uniform weights are applied to vertexes and edges. In Section 3.3, weights coming from the value of the local τ -test will be employed. Graph partitioning algorithms aim to provide a partition that minimizes both the unbalance and the cut of the partition. In our case, the unbalance is not critical, whereas minimizing the cut is essential. The cut of the partition is given by:

cut (P Na ) = k<p cut (C k , C p ) = k<p Ω r ∈C k Ω s ∈Cp ω e (Ω r , Ω s ) (12)
where ω e (Ω r , Ω s ) is the weight of the edge between Ω r and Ω s . The cut is directly linked to the spread of the preconditioned residual Z i , which controls the performance of the algorithm (see Section 2.3 in [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF]) since it drives the cost of the computation of F W i which is operated as:

F W i = F Z i -(F AG)(G AG) -1 G Z i - i-1 j=0 Q j ∆ † j Φ i,j
The matrix (F AG) is computed during the initialization and the last critical point is the computation of:

F Z i = N d s=1 F s Z i
Computing all F s Z i means parallel resolution of local problems (with direct solvers) with multiple righthand sides. The matrix Z i is sparse and the cut of the partition directly defines the number of right-hand sides. Indeed, due to the aggregation, F s needs only to be applied on the columns of Z i which correspond to aggregates containing either Ω s or one of its neighbors Ω s . That is why, the cut of the partition has to be minimized. Aggregation is thus a simple tool to control the cost of the multipreconditioning. The remaining question is, for a given number of aggregates N a , how to choose the partition that leads to the fastest decrease of the error while preserving a low cut. Giving a general answer to this question seems to be out of reach, anyhow it is possible to derive efficient heuristics.

In the simple case of a soft matrix with localized stiff inclusions, it has been observed that the subdomains in the vicinity of the heterogeneity contain more numerical information [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF]. This information was properly caught by the local adaptive test since the search direction proposed by these subdomains were more often selected.

Starting from this observation, two different strategies to adapt the aggregates over the iterations are presented in the following. Both are based on the value of the local adaptive test. In Section 3.2, Section 3.3 and in all numerical tests, the number of aggregates N a is fixed over the iterations. This is not a limitation of the method, but a technical choice since a more efficient memory reuse is obtained in the code by keeping N a constant.

Discrete aggregate adaptation

At each iteration, the N a aggregates are separated into three disjoint groups depending on the value of the adaptive test.

P Na = S c J c U
The set S c := C ∈ P Na | t C i τ contains aggregates that are candidates to be split since, according to the value of the local τ -test, they contain valuable information to achieve convergence. Conversely, the set J c := C ∈ P Na | t C i τ gathers the aggregates that are candidate to be joined. Aggregates that are not in the two previous sets remain unchanged U := P Na \ (J c ∪ S c ). The "much greater" and "much less" relations are defined by two heuristics linked with the cardinal of the aggregate. The proposition

t C i τ is considered to be true if t C i card C ≤ τ . Analogously, the proposition t C i τ is considered to be true if t C i / card C ≥ τ .
Once the aggregates have been classified, various strategies can be used to actually perform the merger and the split.

Join strategies Two strategies are described in Algorithm 2 and 3. The former simply proposes to join aggregates two-by-two if they are connected, i.e. if they contain, at least, two neighboring subdomains. This strategy makes aggregates smoothly evolve and it preserves their connectivity. The latter joins all the connected components of the subgraph associated with the aggregates to be joined. It induces a more massive adaptation while still ensuring the connectivity.

For both algorithms, aggregates are sorted with respect to the value of (t C i / card (C)) C to favor the largest values of the τ -test : aggregates associated with large (t C i / card (C)) are joined first. After the joining step, the number of aggregates has decreased and some free space is available for the split step.

Algorithm 2: Join strategy based on a two-by-two merging.

Input : Set of aggregates selected for join

J i c , criteria (t C i / card (C)) C Output: Set of actually joined aggregates AJ i c 1 Sort J i c in descending order with respect to (t C i / card (C)) 2 Build the subgraph G J = (J i c , Υ C ) 3 AJ i c = {} 4 for C p ∈ J i c do 5 for C q ∈ J i c \ {C p } do 6 if C p and C q are connected then 7 Insert C p ∪ C q in AJ i c 8
Remove C p and C q from J i Input : Set of aggregates selected for join

J i c , criteria (t C i / card (C)) C Output: Set of actually joined aggregates AJ i c 1 Sort J i c in descending order with respect to (t C i / card (C)) 2 Build the subgraph G J = (J i c , Υ C ) 3 Compute the connected components of G J (via Depth First Traversal) 4 AJ i+1 c = {} 5 foreach connected component P of G J do 6 Insert (∪ C∈P C) in AJ i c 7 Remove (C) C∈P from J i c 8 Return AJ i c
Split strategies In the same way, two split strategies are described in Algorithm 4 and 5. Algorithm 4 proposes to bisect aggregates as long as there is free space available (keep in mind that total number of aggregates is fixed). Algorithm 5 enables to split selected aggregates into more than two parts. Algorithm 6 summarizes the adaptation process. For computational efficiency there is a one-iteration delay between the classification and the aggregates update. Indeed, with this delay the local contribution to the preconditioned residual Z p i does not need to be recomputed. The discrete adaptation has the advantage to present a negligible computational cost. However, the evolution of the aggregates may be really slow, even stopped, when the number of aggregates is much smaller than the number of subdomains (see Section 4). This phenomenon has motivated the development of a more continuous process presented hereafter.

Aggregate adaptation via weighted graph partitioning

At each iteration, a new partition (C k ) 1≤k≤Na is computed by a graph partitioning library. The value of the τ -test is used to define the vertexes' weight in the graph, leading to a balanced partition. More precisely, the weights are proportional to 1/t C i in order to correctly adapt aggregates. With this choice, the vertex weight is really high when the local τ -test has detected some "meaningful information". It thus produces aggregates with few subdomains in these areas. In the case of heterogeneous plates with localized inclusions (see Section 4.2), it leads to small aggregates near the heterogeneity for instance. All the subdomains of a given aggregate C are assigned the same vertex weight:

∀Ω s ∈ C, ω v (Ω s ) = 1 1 + t C i Na k=1 t k i ( 13
)
As can be seen in ( 13), the value of the τ -test t C i is normalized with respect to the values coming from all aggregates. When the number of aggregates is small, this strategy is expected to give better results than the discrete one since the adaptation is much less constrained. The computational cost remains small since the number of edges of the subdomain connectivity graph remains low because this graph only concerns short range interactions.

Numerical results

In this section, the performance of the AMPFETI variants are compared. Section 4.1 compares the results obtained by AMPFETIG and AMPFETIL with various numbers of fixed aggregates. A weak scalability study on a checkerboard heterogeneous cube emphasizes that grouping search directions into aggregates provides a better control of the memory footprint without significantly degrading the time to solution. In Section 4.2, the benefit of the two-level adaptive strategy of AMPFETIL is demonstrated with a weak scalability test with localized inclusions. Let us note that FETI and MPFETI algorithms are not taken into account here. A detailed comparison of FETI, MPFETI, AMPFETI showing the good performance of AMPFETI can be found in [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF].

Remarks on the implementation and hardware computing resources

The solvers have been implemented in the finite element suite Z-Set 9.01 . All the computations presented in this section have been performed using the Part 1 of the Occigen supercomputer. This supercomputer is made of Haswell bi-processors (E5-2690V3@2.6 GHz) which are interconnected with an Infiniband FDR 56 Gb/s network. In all configurations, MUMPS solver (version 5.1.2) [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF] is used in association with the BLAS library provided by Intel 17.2 MKL for local solves. The Pardiso direct solver is used to solve the coarse problem. Eigen library2 is used for dense linear algebra. Communication are handled by the OpenMPI library.

Weak scalability study on a highly heterogeneous checkerboard cube

This subsection is dedicated to the study of AMPFETI's scalability. The test-case is a cube with checkerboard heterogeneity pattern. Automatic graph partitioning is employed, which results in numerical difficulties spread all over the structure, making it uninteresting to try adapting aggregates. The performance of AMPFETIG is also mentioned since this variant of the solver seems well suited to that configuration.

For n c ∈ {5, 6, 7, 8, 9, 10, 11}, we consider a set of three-dimensional heterogeneous cubes made of n 3 c identical sub-cubes (see Figure 3). Each sub-cube is discretized with the same ruled mesh made of 110592 eight-node brick elements (c3d8), leading to a total number of approximately n 3 c × 352, 947 degrees of freedom. With this setup, the H/h ratio equals to 48 where h is the diameter of the finite elements and H that of a typical subdomain. The cube is clamped on one face and subjected to a prescribed unitary displacement in the three space directions on the opposite face, all other faces being traction-free. The material behavior is isotropic linear elastic, with a Poisson's coefficient of 0.3 and two values of Young's modulus assigned following a checkerboard pattern in order to obtain a coefficient jump E 1 /E 2 between two adjacent sub-cubes. Two ratios of Young's modulus are used: 10 3 and 10 6 . Finally, an unstructured decomposition in N d = n 3 c subdomains is obtained with a graph partitioning software which leads to interfaces not aligned with the heterogeneity. For a given number of subdomains, the partitioning is computed once and reused for all solvers configurations and for both coefficient jumps. The choice N d = n 3 c , combined with the use of an automatic graph partitioning software leads to a lot of traversing heterogeneities that are known to strongly deteriorate the convergence of the FETI method. Such a configuration is represented in Figure 3 for n c = 6. Six cores are allocated to each subdomain, a shared memory parallelism is used at several steps including (but not limited to) local operators and coarse problem factorization. The study starts from 125 subdomains and goes up to 1728 subdomains which corresponds to a total number of 10368 cores and 576 millions unknowns.

All algorithms are equipped with the best state of the art preconditioner and projector, they use a combination of Primal Schur complement and stiffness scaling. The convergence is triggered when r i / r 0 ≤ . The convergence criterion is set to 10 -6 and the τ -test threshold is set to 10 -2 . The influence of the choice of the preconditioner, projector and threshold is discussed at length in [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF]. Several numbers of aggregates have been tested N a ∈ {32, 64, 125}. Table 1 summarizes the different configurations.

The results are shown in Figures 45. The letter G (resp. L) stands for AMPFETI with the global τ -test (resp. local τ -test), the appended number indicates the number of aggregates. The total wall time includes the setup, the factorization of local operators and the convergence loop of the algorithm. For both AMPFETIG and AMPFETIL, and whatever the Young's modulus ratio, the same behavior is observed.

• The number of iterations needed to achieve convergence slightly increases with the number of subdomains. This behavior is not uncommon for such high level of heterogeneity and automatic decomposition. In GENEO's approaches, it is often observed that the number of modes to be eliminated per subdomains increases with the number of subdomains. In AMPFETI approaches, the bad modes must be captured through extra iterations, at least one iteration per bad mode, before superconvergence can be observed.

• Regarding the time to solution, a change of slope is clearly visible after 729 subdomains (n c = 9). This change of slope is induced by two phenomena. First, the time spent in building and solving the coarse problem strongly increases when the number of subdomains becomes large in our implementation. For now, G AG is communicated to all "subdomains" and all of them perform the factorization with a multithreaded solver (redondant work). A global communication with a significant amount of data happens here. The implementation of a better strategy, inspired from [START_REF] Hapla | TFETI Coarse Space Projectors Parallelization Strategies[END_REF], will be considered in the future. Also, a superlinear increase of the cost of all MPI_Allreduce operations involved by the full reorthogonalization is observed (even if the increase of the search space size is limited).

• As expected, the size of the search space decreases with the use of smaller number of aggregates.

For the largest test case (N d = 1728) and whatever the heterogeneity, the search space generated with 32 aggregates is approximately 15% smaller than the one produced with 125 aggregates. The increase of the time to solution is limited since the additional time induced by the higher number of iterations is somehow compensated by time saved in the orthogonalization.

• The increase of the number of iterations with respect to N d is slightly faster when the Young's modulus ratio is 10 6 . It was expected because n c also increases leading to more and more material coefficent jumps in the global problem and thus deteriorating the condition number.

Weak scalability study on a three dimensional elastic plate with stiff inclusions

While the previous section has shown the benefit of gathering search directions into aggregates, this section highlights the advantage of the two-level adaptation. Only AMPFETI with the local τ -test is shown here since it is the only algorithm that is concerned with this adaptation. The pattern is defined by a three-dimensional heterogeneous plate made of a soft matrix and stiff inclusions. In this pattern, two hundred inclusions are inserted in two strips 0 ≤ x ≤ 0.2 x max and 0.6 x max ≤ x ≤ 0.8 x max (see Figure 6). The matrix behavior is isotropic elastic linear with a Poisson's coefficient of 0.45. All spherical inclusions have a Poisson's coefficient of 0.3. The Young's modulus of the inclusions are randomly selected using uniform distribution, such that the maximum Young modulus ratio E inclusion /E matrix is 10 6 . The plate is clamped on one side and subjected to a prescribed unitary displacement in the three space directions on the opposite side, all other faces remain traction-free. The pattern is discretized with a ruled mesh made of 125 × 125 × 10 eight-node brick elements (c3d8), leading to a number of 174, 636 nodes. To perform the weak scalability study, this pattern is repeated in the directions x and y, and an unstructured decomposition is obtained with a partitioning software. The number of subdomains per pattern is 120. Thus, the study starts from 2 × 2 patterns (480 subdomains) and goes up to 7×7 patterns (5880 subdomains). Since subdomains are really small in this test case (about 5, 000 unknowns per subdomain), only one core is attributed to each one. Table 2 summarizes the different configurations. All algorithms are equipped with stiffness scaling and Dirichlet preconditioner for S and A. The convergence is triggered when r i / r 0 ≤ . The convergence criterion is set to 10 -6 and the τ -test threshold is set to 10 -2 . Several numbers of aggregates have been used N a ∈ {32, 64, 125, 216}. For each number of aggregates, the adaptation based on weighted graph partitioning and two combinations of discrete split and join strategies have been tested. Table 3 gives the correspondence between the legend and the adaptation strategy presented in Section 3. The letter L stands for the non-adaptive AMPFETI with local test.

The results are summarized in Tables 4567. Figure 8 compares the scalability of the L 2 variant whereas Figure 7 presents a focus on the 5880-subdomain case.

We observe that for all cases the convergence of the two-level adaptive variants is faster than the fixed aggregates variant (L). As shown in Tables 4 andTables 5, when the number of aggregates is much smaller than the number of subdomains, all discrete adaptation variants remain inactive and do not provide any improvement. This is not the case for the continuous variant. As seen in Figure 8, the L 2 variant always leads to fewer iterations than the fixed aggregate configuration. For the case with 5880 subdomains and 64 aggregates, the gain is about 22%. More, the search space size for both configurations is approximately equal. Because of the small size of the subdomains (about 5000 unkowns), the gain in terms of number of iterations has little effect on the time to solution. Indeed, the time spent in local solves is much less than the time spent in the coarse problem. So even if a gain is obtained in terms of iteration, it is not so much visible in the computational time. This setup was chosen to be able to reach a large number of subdomains without consuming too much computational resources. The search space size reduction is about 10% between N a = 32 and N a = 216. The aggregate process leads to moderate gain here because, contrary to the checkerboard test case, the local τ -test of AMPFETIL does not saturate here but remains "frugal".

N d #DOFs
As presented in Figure 7, with the two-level adaptation method, the algorithm usually keeps selecting search directions a little longer, which means more relevant search directions are generated. aggregates. In Figure 9b-9f, each dot represents a subdomain. These dots are situated at the gravity center of subdomains and they are linked if they belong to the same aggregate. Figure 9a is an above view of the geometry allowing to compare the aggregates evolution with the position of the stiff inclusions.

Starting from a uniform distribution of aggregates, as soon as the second iteration, the algorithm generates "small" aggregates close to the stiff inclusions. After iteration 11, the distribution obtained is the final one. In Figure 9f, one can observe that close to the inclusions, allmost all aggregates only contain one subdomain. Conversely, aggregates localized in the homogeneous part of the problem gather several, and sometimes many, subdomains.

Applications to real engineering problems

This section presents some typical real engineering problems where the use of the AMPFETI solver is essential to achieve convergence.

Numerical homogenization of the mechanical behavior of solid propellants

Solid propellants are energetic materials composed of an organic matrix and numerous metallic inclusions. Due to the large dispersion of the sizes of the particles, this type of material brings several space scales into play. The simulation of a Representative Volume Element (RVE) with the finite element method leads to many unknowns, making the use of iterative solvers mandatory. However, the numerical homogenization of this type of material with domain decomposition methods is a challenging problem since it gathers almost all classical hard points that penalize the convergence of iterative solvers. First, the material is highly heterogeneous, the jump in Young modulus between the matrix and the metallic inclusions is approximately 10 5 . Because of the high density of inclusions in the RVE, heterogeneity is very frequently misplaced with respect to the domain decomposition interface. Also, elements of the mesh of the matrix that are located between two close inclusions often exhibit poor quality factors which degrade even more the condition number of the linear system to be solved. More, the organic matrix is almost incompressible, and mixed pressure-displacement-volume variation finite element needs to be used [START_REF] Al Akhrass | Integrating a logarithmicstrain based hyperelastic formulation into a three-field mixed finite element formulation to deal with incompressibility in finite-strain elastoplasticity[END_REF] (if a linear elastic behavior is assumed for the matrix, the Poisson coefficient is 0.499). This is a typical case where the AMPFETI solver is needed to achieve convergence. As an example, we present some results on the linearized problem. Figure 10a shows the finite element mesh and the domain decomposition for one distribution of particles. The finite element mesh contains 5,494,528 quadratic tetrahedrons and the global system has 22 millions of degrees of freedom; 448 subdomains are used for the partitioning. The computation is run on 448 cores (MPI parallelism only). Figure 11 shows the convergence of the solver and the generated search space for various configurations. The curves G448 and L448 correspond to the original AMPFETI (with one subdomain per aggregate, G means global τ -test, L means local τ -test). Although the inclusions are spread out over all the RVE, all local variants generate a smaller search space than their global counterpart. It is especially visible when comparing G448 with L448. The use of aggregates allow to reduce the search space size without degrading too much the convergence. Here, the second level of adaptation (L 2 112 and L 2 56) leads to the same results as the one level (L112 and L56). It may be explained by the fact that inclusions are spread out over all the RVE and not localized.

High fidelity simulation of a multiperforated aircraft combustion chamber

Aircraft combustion chambers are complex systems subjected to strong thermomechanical loading. To estimate the fatigue lifespan, several full load cycles need to be simulated. Thanks to AMPFETI, the simulation can be performed in a reasonable CPU time, while taking into account the anisothermal elastoviscoplastic highly nonlinear material and the fine description of the geometry. The material used is modeled with a complex elasto-viscoplastic law. This law has two potentials, the first potential is elastoviscoplastic, it includes a Norton flow, two hardening mechanisms (one kinematic and one isotropic).

The second potential is viscoplastic with a Norton flow. All material parameters are temperature dependant. Also, a complete set of 1200 micro-perforations is modeled and their influence on the stress and temperature field can be analyzed. Figure 12a shows the domain decomposition. A zoom on the finite element mesh close to the perforations is visible in Figure 12b. The finite element mesh contains quadratic tetrahedrons and hexahedrons, the global system has approximately 25 million degrees of freedom. There are 576 subdomains, the computation is run on 3456 cores (hybrid MPI-OpenMP parallelism). Because of the material nonlinearity, a Newton-Raphson scheme is needed and AMPFETI is used as a linear solver for the successive tangent linear systems. Each Newton iteration lasts approximately one minute, this time includes the resolution of the tangent system and the integration of the nonlinear material law. As an example, the Von mises stress field at some instant of the cycle is shown in12c.

Figure 13 shows the convergence of the solver and the generated search space for various configurations for the first tangent system. The curve L576 corresponds to the original AMPFETI (with one subdomain per aggregate). Again, the use of aggregates allows a significant reduction of the search space size without degrading too much the convergence. Also, the second level of adaptivity (L 2 288 and L 2 144) leads approximately to the same results as their one level counterparts.

Conclusion and perspectives

This paper presents strategies to reduce the memory cost of the Adaptive Multipreconditioned FETI method.

First, the gathering of search directions into aggregates allows a better control of the memory requirements of the solver since it limits the search space growth per iteration. It is especially useful for large scale ill conditioned systems.

Then, a second level of adaptation has been introduced to improve the convergence by updating aggregates at each iteration. It is especially efficient when the computation complexity is localized in few subdomains and when the number of multipreconditioning aggregates is small in comparison to the number of subdomains. The update of the multipreconditioning aggregates makes use of the local τ -test. Two evolution strategies have been derived. The first one is based on aggregates classification into three sets (split candidates, join candidates, unchanged); the update of aggregates looks like a discrete system with a limited number of states. This strategy has the advantage to be costless, but it may be inadequate when the number of aggregates is too small. The second one, based on weighted graph partitioning, is more flexible but has a higher computational cost.

The benefits of these modifications have been assessed with two weak scalability studies going up to 10368 cores and 576 millions of unknowns. Finally, two engineering simulations leading to ill-conditioned systems have been presented. In these simulations, the use of AMPFETI was essential to achieve convergence.

One perspective of this work is to replace the sum in the gathering process by a weighted sum and to look for optimal coefficients. Other perspectives are to associate multipreconditioning with classical deflation and restarting techniques and to propose adaptive strategy for the non symmetric case.

  and use the definition[START_REF] Bovet | Two-level adaptation for adaptive multipreconditioned FETI[END_REF] of Z i+1 at line 13.

Figure 1 :

 1 Figure 1: Cubic mesh with regular decomposition

Figure 2 :

 2 Figure 2: Illustration of a subdomain's connectivity graph. The graph is decomposed in two colors (blue and red). The dashed line represents the separation. The cut of the partition is ω e (Ω 1 , Ω 3 ) + ω e (Ω 2 , Ω 3 ).

c 9 break 10 Algorithm 3 :

 9103 Return AJ i c Join strategy based on the connected components of the subgraph.

Algorithm 6 : 4 5 Local test with the newly formed aggregates ( 12) 6 Classify, U i+1 ) 7 AMPFETI

 64567 AMPFETI with discrete aggregate adaptation based on the local τ -test. Input: Initial aggregate partition split and join strategies.1 AMPFETI initialization ( 1 -6) 2 while r i > do 3 AMPFETI beginning of iteration( 8 -11 ) Adapt aggregates according to (S i c , J i c , U i ) end of iteration( 13 -14 )

Algorithm 7 : 4 Call the graph partitioner and adapt the aggregates 5 6

 7456 Aggregate adaptation via weighted graph partitioning based on the local τ -test. Input: Initial aggregate partition 1 AMPFETI initialization ( 1 -6) 2 while r i > do 3 AMPFETI start of iteration ( 8 -11 )Local test with the newly formed aggregates[START_REF] Farhat | Extending substructure based iterative solvers to multiple load and repeated analyses[END_REF] Compute the vertex weights (ω(Ω s )) 1 s N d 7 AMPFETI end of iteration( 13 -14 )

  (a) Checkerboard cube. Red and blue areas correspond to the two different materials. (b) Automatic decomposition. Each color represents a different subdomain.

Figure 3 :

 3 Figure 3: Heterogeneous cube (configuration with N d = 216, n c = 6).

  (c) Minimization space size.

Figure 4 :

 4 Figure 4: cube, weak parallel scalability (E 1 /E 2 = 10 3 ): wall time, number of iterations and minimization space size.

  (c) Minimization space size.

Figure 5 :

 5 Figure 5: Checkerboard cube, weak parallel scalability (E 1 /E 2 = 10 6 ): wall time, number of iterations and minimization space size.

  (a) Heterogeneous 3D plate pattern. (b) Test case with 2 × 2 repetition.

Figure 6 :

 6 Figure 6: Weak scalability with localized inclusions.

  Convergence results with 125 aggregates.

  (d) Convergence results with 216 aggregates.

Figure 7 :

 7 Figure 7: Weak scalability test with localized inclusions, focus on the results for the case with 5880 subdomains.

32 (

 32 c) Minimization space size.

Figure 8 :

 8 Figure 8: Weak scalability test with localized inclusions, focus on the results of the graph based update of aggregates: wall time, number of iterations and minimization space size.

  (a) Above view of the problem.

  Aggregates at iteration 6.

  Aggregates at iteration 11 (final).

Figure 9 :

 9 Figure 9: Aggregates evolution for variant B with 480 subdomains and 120 aggregates. Each dot represents a subdomain. Dots are linked if they belong to the same aggregate.

  (a) Finite element mesh of the RVE, the red and yellow colors correspond to two different types of inclusions. (b) Domain decomposition in 448 subdomains. (c) Von Mises stress field, the white color characterises low stresses area.

Figure 10 :

 10 Figure 10: Numerical homogenization of the mechanical behavior of solid propellants.

  Evolution of the total search space size.

Figure 11 :

 11 Figure 11: Solid propellant: convergence and search space size. The local preconditioner is the Dirichlet one with a k-scaling. The user parameter τ is 0.01.

  (a) 576-subdomain decomposition. (b) Finite element mesh. (c) Von mises stress fields

Figure 12 :

 12 Figure 12: Multiperforated aircraft combustion chamber simulation.

  (b) Evolution of the total search space size.

Figure 13 :

 13 Figure 13: Multiperforated aircraf combustion chamber: convergence and search space size for the first tangent system. The local preconditioner is the Dirichlet one with a k-scaling. The user parameter τ is 0.1.

  • A global adaptive test (called Global τ -test) in which either all search directions are kept (leading to a full MPCG iteration), or they are all reduced to one direction (leading to a classical PCG). • A local adaptive test (called Local τ -test) in which each search direction is considered individually. This variant allows a much finer selection of search direction candidates. But the local test is significantly more expensive to compute than the global test.

Algorithm 4 :

 4 Split strategy based on aggregates bisection. Split strategy based on a uniform partitioning of aggregates.

		Input : Set of aggregates selected for split S i c , criteria (t C i card (C)) C , free space Output: Set of actually split aggregates AS i c
	1 Sort S i c in ascending order with respect to (t C i card (C))
	2 AS i+1 c 3 for C p ∈ S i = {} c do
	4	if freespace is available then
	5	Append bisect(C p , 2) to AS i c
	6	Update the freespace
	7 Return AS i c
	Algorithm 5: Input : Set of aggregates selected for split S i c , criteria (t C i card (C)) C , F (free space) Output: Set of actually split aggregates AS i c
	1 Sort S i c in ascending order with respect to (t C i card (C))
	2 AS i c = {} 3 P = ones(card(S i c ))
	4 while F > 0 (free space is available) do
	5	for C p ∈ S i c do

6 if P (p) < card(C p ) then 7 P (p) ← P (p) + 1 8 F ← F -1 (update the free space) 9 for C p ∈ S i c do 10 Split C p into P (p) parts 11 Append the parts to AS i c 12 Return AS i c

Table 1 :

 1 Checkerboard cube, weak parallel scalability: configurations.

	nc	N d	#DOFs total #DOFs on interface #cores
	5	125	41,992,563	2,645,043	750
	6	216	72,412,707	4,785,987	1296
	7	343	114,818,259	7,752,734	2058
	8	512	171,199,875	11,641,947	3072
	9	729	243,548,211	17,084,832	4374
	10 1000	333,853,923	23,613,033	6000
	11 1331	444,107,667	32,187,410	7986
	12 1728	576,300,099	42,008,133	10368

Table 2 :

 2 Three dimensional plate with stiff inclusions, weak parallel scalability: studied configurations.

		total #cores
	480	2,079,033	480
	1080	4,665,408	1080
	1920	8,283,033	1920
	3000	12,931,908	3000
	4320	18,612,033	4320
	5880	25,323,408	5880

Table 3 :

 3 Three dimensional plate with stiff inclusions, weak parallel scalability: tested combinations.

Table 4 :

 4 Weak scalability with heterogeneous inclusions: summary of results for the case with 32 aggregates. The user parameter τ is set to 0.01.

	N d	480	1080 1920 3000 4320	5880
			Number of iterations
	L	83	118	184	267	336	389
	A	83	118	184	267	336	389
	B	83	118	184	267	336	389
	L 2	80	96	155	217	274	342
			Search space dimension
	L	1014 2076 3571 5603 7689 10381
	A	1014 2076 3571 5603 7689 10381
	B	1014 2076 3571 5603 7689 10381
	L 2	1081 1995 3520 5607 7712 10548
			Time to solution (s)
	L	8.865 33.9 191.0	249	455	911
	A	8.777 30.0 132.1	207	445	945
	B	8.79	29.5 200.5	208	463	922
	L 2 9.013 28.7	81.9	201	440	843
	N d	480	1080 1920 3000 4320 5880
			Number of iterations
	L	68	79	124	161	190	259
	A	65	79	124	161	190	259
	B	65	79	124	161	190	259
	L 2	64	69	106	140	168	202
			Search space dimension
	L	1047 2130 3631 5732 7821 10723
	A	1062 2130 3631 5732 7821 10723
	B	1092 2130 3631 5732 7821 10723
	L 2 1206 2204 3645 5682 7774 10722
			Time to solution (s)
	L	10.08 30.26 98.5	208	444	791
	A 9.075 28.97 85.1	213	390	790
	B 9.072 33.56 85.6	212	386	756
	L 2 10.09 28.98 85.0	211	444	750

Table 5 :

 5 Weak scalability with heterogeneous inclusions: summary of results for the case with 64 aggregates. The user parameter τ is set to 0.01.

	N d	480	1080	1920 3000 4320	5880
			Number of iterations
	L	56	59	87	108	126	153
	A	52	55	87	108	126	153
	B	51	56	87	108	126	153
	L 2	56	53	81	102	114	130
			Search space dimension
	L	1229	2339	3732 6188 8003 11065
	A	1400	2222	3732 6188 8003 11065
	B	1361	2341	3732 6188 8003 11065
	L 2	1271	2443	4038 5927 7939 10781
			Time to solution (s)
	L	10.33 30.61 105.1	201	413	741
	A	10.85 27.02	79.5	192	361	750
	B	10.49	30.3	80.9	234	432	746
	L 2 10.44 30.12	97.7	191	418	757

A focus on aggregates evolution In order to "feel" how the algorithm behaves, we propose to visualize the aggregates evolution given by the variant B on the test case with 480 subdomains with 120

Table 6 :

 6 Weak scalability with heterogeneous inclusions: summary of results for the case with 125 aggregates. The user parameter τ is set to 0.01.

	N d	480	1080 1920 3000 4320	5880
			Number of iterations
	L	51	47	72	86	109	117
	A	48	42	66	88	102	117
	B	47	41	71	86	102	117
	L 2	49	46	65	83	89	96
			Search space dimension
	L	1502	2697 4229 6621 8809 11252
	A	1449	2599 4181 6367 8761 11252
	B	1522	2671 4021 6412 8761 11252
	L 2	1529	2800 4313 6374 8413 11118
			Time to solution (s)
	L	11.66 33.86 85.4	199	427	756
	A	11.04 29.52 83.5	191	476	768
	B	11.51	31.4	78.5	198	418	776
	L 2 11.51	31.1	83.9	288	419	738

Table 7 :

 7 Weak scalability with heterogeneous inclusions: summary of results for the case with 216 aggregates. The user parameter τ is set to 0.01.

http://www.zset-software.com/

http://eigen.tuxfamily.org/

Acknowledgements

This work benefited from the support of the project SEMAFOR ANR-14-CE07-0037 of the French National Research Agency (ANR). This work was granted access to the HPC resources of CINES under the allocation 2017-A0040910442 made by GENCI. The author would also like to thanks V. Chiaruttini for providing the engineering test cases presented in Section 5.