Thiago Milanetto

Schlittler • Régis

Fully scalable implementation of a volume coupling scheme for the modeling of multiscale materials

Keywords: Multimodel, Multiscale, Arlequin method, FETI method, Mesh intersection search

We present in this paper a new implementation of a multi-scale, multi-model coupling algorithm, with a proposed parallelization scheme for the construction of the coupling terms between the models. This allows one to study such problems with a fully scalable algorithm on large computer clusters, even when the models and/or the coupling have a high number of degrees of freedom. As an application example, we will consider a system composed by an homogeneous, macroscopic elasto-plastic model and an anisotropic polycrystalline material model, with a volume coupling based on the Arlequin framework.

Introduction

Coupling methods, in general, allow the study of multiscale systems, taking into account the different physical processes of each scale -for example, the coupling between continuous and granular mediums [START_REF] Wellmann | A two-scale model of granular materials[END_REF], and of heterogeneous wave propagation solvers [START_REF] Marchais | Representation of localized phenomena in dynamics using multi-scale coupling[END_REF][START_REF] Ghanem | Arlequin framework for multimodel, multi-time scale and heterogeneous time integrators for structural transient dynamics[END_REF]. These methods can be classified in several different types, depending on the desired application. Methods like the VMS [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods[END_REF] and the HMM [START_REF] Engquist | The Heterognous Multiscale Methods[END_REF], and similar others, are used when the parameters of the macro-scale model are not known, but depend on the micro-scale over the whole domain. They can be classified as embedding methods, with the VMS enriching the macro-scale with the micro-scale model over element patches, and the HMM altering the quadrature evaluation of the weak formulation. If, on the other hand, the macro-scale quantity of interest depends on the micro-scale only over a part of the former's domain, more local methods are used. Examples include the non-overlapping domain decomposition methods [START_REF] Alart | A nonlinear Domain Decomposition formulation with application to granular dynamics[END_REF]. They are derived from domain decomposition techniques developed to solve numerically PDEs over large-scale computer clusters [START_REF] Dolean | An Introduction to Domain Decomposition Methods[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF], and they are used when the coupling is done over an interface between the models. These models use the following formulation: find (u 1 , u 2 , Φ) ∈ W 1 ×W 2 ×M such that

a 1 (u 1 , v 1 ) + c (Φ, v 1 ) = 1 (v 1 ) , ∀v 1 ∈ W 1 ; (1a) a 2 (u 2 , v 2 ) -c (Φ, v 2 ) = 2 (v 2 ) , ∀v 2 ∈ W 2 ; (1b) c (ψ, Π 1 u 1 -Π 2 u 2 ) = 0, ∀ψ ∈ M, (1c)
where the operators Π l : W l → M project functions from W l into their restriction on the interface. The weak formulations of the coupled models (first and second equations) are altered by adding a coupling operator c(•, •) and a Lagrange multiplier Φ, defined over a mediator space M and associated to the interface between the models. Let us note as W l the spaces associated to each model l. The last equation guarantees that the solutions of the systems, when projected into M, are continuous under the coupling operation. The formulation Eq. ( 1) can be solved iteratively by, for example, enforcing strongly the continuity condition, or by computing the Lagrange multiplier as an intermediary step [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF].

For overlapping domains, still in the context of a limited domain coupling, methods based on volume couplings can be used, such as the Arlequin framework [START_REF] Ben Dhia | Problèmes mécaniques multi-échelles: la méthode Arlequin (written in French)[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF]28] and the bridging domain method [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF]. They keep the same formulation as Eq. ( 1) but differ from the interface methods by defining the term c(•, •) as a volume coupling over the overlapping domains. These methods differ mainly on the choice of this coupling term. In both cases, c(•, •) is built using an intermediary mesh, and a specific effort must be set on integrating functions borne by incompatible meshes: that related to each of the models, and that of the coupling operator c(•, •).

The problem of meshing geometrical intersections is also present in methods where a volume is cut by an interface, such as with Nitsche's method [START_REF] Fernández-Méndez | Imposing essential boundary conditions in mesh-free methods[END_REF][START_REF] Juntunen | Nitsche's mthod for general boundary conditions[END_REF][START_REF] Embar | Imposing Dirichlet boundary conditions with Nitsche's method and splinebased finite elements[END_REF] or the eXtended Finite Element Method [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Chessa | An extended finite element method for two-phase fluids: flow simulation and modeling[END_REF][START_REF] Díez | A stable extended FEM formulation for multi-phase problems enforcing the accuracy of the fluxes through Lagrange multipliers[END_REF]. In the former method, two overlapping meshes are considered: one representing an object, and another representing its surroundings -which is altered to exclude the overlapping between both meshes. To do this operation, one needs to identify the geometric intersections between the interface of the object mesh and the background mesh, and re-mesh the intersecting elements of the latter. In the latter method, functions with discontinuities (or with discontinuous derivatives) along a given interface have to be integrated over a volume mesh. Both these methods and the overlapping coupling techniques therefore share similar characteristic features in terms of integration, which are discussed in more detail in [START_REF] Massing | Efficient implementation of finite element methods on non-matching and overlapping meshes in 3d[END_REF].

Generally speaking, the coupling step is not parallelized for methods following the formulation Eq. ( 1), since the mediator space is, in many cases, small when compared to the models' spaces associated to it. By consequence, while the steps associated with each model scale well when parallelized, the serial coupling step breaks this scalability of the algorithm as a whole. Previous works [START_REF] Ben Dhia | Multimodeling of multi-alterated structures in the Arlequin framework[END_REF][START_REF] Elkhodja | Approches de structures complexes dans des cadres adaptés de la méthode Arlequin[END_REF] using the coupling algorithms and parallelism focused on simulating several couplings between a single global model and many local models, with each coupling being associated to one processorhence each coupling is serial. Here, we will use the same theoretical framework as presented in ref. [START_REF] Ben Dhia | Multimodeling of multi-alterated structures in the Arlequin framework[END_REF][START_REF] Elkhodja | Approches de structures complexes dans des cadres adaptés de la méthode Arlequin[END_REF], but we will focus on an implementation with parallelized and scalable coupling step, possible due to a detailed work on all steps of the implementation.

As an example of this implementation, we will study a coupled system formed by a macroscopic and homogeneous elastic model and a microscopic anisotropic elastic model, representing a polycrystalline material. While we will focus on this example here, this framework -and thus this implementation -can be applied to other cases involving multi-scale physics, such as the ones cited above. Similarly, we insist that, while we focus here in the Arlequin framework, the parallelization scheme presented here is applicable to any coupling methods using the formulation Eq. [START_REF]Code Arlequin[END_REF].

This paper is structured on the following manner. We will start by defining our physical problem, and by recalling the Arlequin framework and the FETI domain decomposition method [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF], which was adapted to solve the Arlequin problem in refs. [START_REF] Ben Dhia | Multimodeling of multi-alterated structures in the Arlequin framework[END_REF][START_REF] Elkhodja | Approches de structures complexes dans des cadres adaptés de la méthode Arlequin[END_REF]. We will follow this with a discussion on the coupling operators -including our implementation of a parallel intersection search algorithm, needed to assemble them. This part should be seen as the main contribution of this article, ensuring the actual scalability (in practice) of a known numerical coupling method (the FETI / Arlequin solver). Finally, we will present in the last section our numerical results, including the scalability tests of the algorithms and an application case with a crack test involving homogeneous and heterogeneous elasticity models.

2 Arlequin framework

Formulation

We will consider here two linear elastic models, each associated to bounded regular domains Ω 1 and Ω 2 , and which overlap over a region Ω 12 . The latter is decomposed into two non-overlapping regions,

Ω 12 = Ω c 12 ∪ Ω f
12 , with Ω c 12 defining the coupling domain between the two models (see Fig. 1). Furthermore, let us note as Γ u and Γ F , respectively, the clamped region of Ω 1 , and its region where a field of surface density forces F is applied.

The functional spaces associated to each model are

W 1 = v ∈ H 1 (Ω 1 ) ; v = 0 on Γ u , W 2 = v ∈ H 1 (Ω 2 ) .
The internal and external works a l (u l , v l ) and l (v l ) are weighted by weight parameter functions (α 1 , α 2 ) de-fined in Ω 1 and Ω 2 , respectively. These weights guarantee the partitioning of the energy between the models on the coupling region, and they follow the relations

α l ≥ 0 in Ω l , l = 1, 2, (2)
α l = 1 in Ω l \ Ω 12 , α 1 + α 2 = 1 in Ω 12 .
The weighted works are defined as

∀ (u 1 , v 1 ) ∈ W 1 ×W 1 , (3)
a 1 (u 1 , v 1 ) = Ω1 α 1 σ (u 1 ) : ε (v 1 ) dΩ, 1 (v 1 ) = Ω1 α 1 f • v 1 dΩ + Γ F α 1 F • v 1 dS, and
∀ (u 2 , v 2 ) ∈ W 2 ×W 2 , a 2 (u 2 , v 2 ) = Ω2 α 2 σ (u 2 ) : ε (v 2 ) dΩ, 2 (v 2 ) = Ω2 α 2 f • v 2 dΩ.
In these equations, f denotes the field of volume density forces applied on both model domains. σ (v l ) and ε (v l ) are the linearized strain and stress tensors, associated to the displacement field v l , and linked through Hooke's law.

The coupled problem involving the two systems follows the formulation presented in Eq. ( 1), with the coupling operator being defined over the coupling domain Ω c

12 . The mediator space, noted by M, is defined as

M = Φ ∈ H 1 (Ω c 12
) . The coupling terms c (•, •) link this mediator space to either of the models spaces, or to itself (in the case of the third equation of Eq. ( 1)) and it follows

∀ (Ψ, Φ) ∈ M × V, V ∈ {W 1 , W 2 , M} c (Ψ, Φ) = Ω c 12 κ ε (Ψ) : ε (Φ) + 1 e 2 Ψ • Φ dΩ. ( 4
)
The parameters κ and e are of the order of magnitude of the material's rigidity and the width of the coupling domain. The coupling term is, then, analogous to a stiffness.

The domains Ω 1 and Ω 2 can be discretized by associating to them, respectively, the meshes T 1 and T 2 , each one with n 1 and n 2 degrees of freedom (DoF's). The formulation Eq. (1) can then be written as

K 1 u 1 + C T 1 Φ = F 1 , (5a)
K 2 u 2 -C T 2 Φ = F 2 , (5b)
C 1 u 1 -C 2 u 2 = 0. ( 5c
)
Following standard notation, K l and F l , l = 1, 2, are the matrices and vectors representing the internal and external virtual works of the model l, with dimensions n l × n l and n l , respectively. The coupling terms C l are n m × n l matrices, where n m is the number of DoF's of M. Each element of this matrix is given by

C l,ij = Ω c 12 κ ε (ν m i ) : ε ν l j + 1 e 2 ν m i • ν l j dΩ, (6)
where ν m i and ν l j are form functions associated respectively to the coupling region and the model l.

If the models' and the mediator form functions can be interpolated exactly by a common function space (as we will do here, in Sec. 4), the coupling matrices C l can be written as the products of a n m × n m coupling matrix, C m and n m × n l interpolation matrices P l :

C l = C m P l . (7)
This new coupling matrix is defined only in the mediator space M, and is given by

C m,ij = Ω c 12 κ ε (ν m i ) : ε ν m j + 1 e 2 ν m i • ν m j dΩ, (8)
where now all the form functions are associated to the coupling region.

3 Solving the Arlequin problem: FETI solver

The system (5) can be solved as a monolithic problem, but this approach has the limitation of not allowing the usage of the proper solvers of the super-imposed models. One can think, for example, of the coupling between models such as linear / nonlinear, deterministic / stochastic, continuum / atomistic ... . Some works have been proposed to solve this problem, adapting domain decomposition methods such as the FETI [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Ben Dhia | Multimodeling of multi-alterated structures in the Arlequin framework[END_REF][START_REF] Elkhodja | Approches de structures complexes dans des cadres adaptés de la méthode Arlequin[END_REF] or the LATIN [28] methods. In the context of these adaptations, each domain of the decomposition corresponds to a different model domain, Ω l , and the domain interfaces correspond to the coupling region Ω c 12 . We focus in this article on a Arlequin solver based on the FETI method. This domain decomposition method, essentially, solves a given model inside each of the domains of the decomposition, and then calculates corrections due to the continuity between these domains. It has a formulation similar to the one presented in Eq. ( 1), with the differences that the binary matrices connecting the domains (noted as B i , in ref. [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]) are exchanged by the coupling matrices C i .

An application of this method to the Arlequin problem is presented in refs. [START_REF] Ben Dhia | Multimodeling of multi-alterated structures in the Arlequin framework[END_REF][START_REF] Elkhodja | Approches de structures complexes dans des cadres adaptés de la méthode Arlequin[END_REF], where the different models and the model couplings correspond, respectively, to the domains of the decomposition and their interfaces. We will present in this section a summarized version of this application, restricted to two coupled models. This will be followed by more detailed description of a solver for the coupling corrections, with a strong focus on which of its components affect the solver's scalability.

FETI solver for well conditioned matrices

If both matrices K 1 and K 2 are well conditioned, the system (5) can be rewritten in the following manner. By isolating the displacements in the first and second equations, we have

u 1 = u 0 1 -K -1 1 C T 1 Φ, (9)
u 2 = u 0 2 + K -1 2 C T 2 Φ, (10)
where u 0 1 and u 0 2 are the solutions of the decoupled models:

K l u 0 l = F l , l = 1, 2. (11)
Substituting Eqs. ( 9) and (10) into Eq. (5c), we obtain a system depending only on the Lagrange multiplier Φ:

C 1 K -1 1 C T 1 + C 2 K -1 2 C T 2 • Φ = C 1 u 0 1 -C 2 u 0 2 (12) AΦ = b. (13)
The solution (u 1 , u 2 , Φ) can then be found by first solving the decoupled systems [START_REF] Díez | A stable extended FEM formulation for multi-phase problems enforcing the accuracy of the fluxes through Lagrange multipliers[END_REF], then solving the coupled system [START_REF] Dolean | An Introduction to Domain Decomposition Methods[END_REF], and finally adding the coupling corrections to the models through Eqs. ( 9) and (10).

3.2 FETI solver with a singular matrix, K 2

The method above must be altered if one of the models has a singular matrix K l (for example, in the context of linear elasticity, if the model has no Dirichlet boundary conditions). Let us suppose here that the second model falls in this case and that K 2 is singular. Eq. (5b) still admits solutions if K 2 is consistent. In this situation, Eq. ( 10) is substituted by

u 2 = K + 2 F 2 + C T 2 Φ + R 2 α = u 0,+ 2 + K + 2 C T 2 Φ + R 2 α ( 14
)
where

K + 2 is the pseudo-inverse of K 2 (a n 2 ×n 2 matrix such that K 2 K + 2 K 2 = K 2 ). The matrix R 2 is a n 2 × n RB 2
matrix, whose columns form a basis of K 2 's null space, with n RB vectors. α is a n RB 2 vector defining a linear combination of this basis. In the context of the elastic models that we are considering here, K 2 's null space is formed by the rigid body modes that are left "free" due to the lack of Dirichlet boundary conditions on the model.

The system (5b) admits at least one solution of the form [START_REF] Elkhodja | Approches de structures complexes dans des cadres adaptés de la méthode Arlequin[END_REF] if and only if the term F 2 + C T 2 Φ is orthogonal to K 2 's null space. This condition is represented by the relation

R T 2 F 2 + C T 2 Φ = 0. ( 15
)
Substituting Eq. ( 9) and ( 14) into Eq. (5c), we obtain a system depending on the Lagrange multiplier and the null space coefficients α:

A + Φ + R I 2 α = b + , (16)
where

A + = C 1 K -1 1 C T 1 + C 2 K + 2 C T 2 , (17)
R I 2 = C 2 R 2 , and (18)
b + = C 1 u 0 1 -C 2 u 0,+ 2 . ( 19
)
Together with [START_REF] Embar | Imposing Dirichlet boundary conditions with Nitsche's method and splinebased finite elements[END_REF], this equation can be written in matrix form as

A + R I 2 R I 2 T 0 • Φ α = b + -R T 2 F 2 . ( 20
)

Projected conjugate gradient algorithm

The system (20) is symmetric and non-singular, and hence it admits an unique solution (Φ, α). Still, this system is also indefinite, so an iterative method such as the conjugate gradient algorithm (CG) cannot be directly applied. Instead, a projected CG algorithm can be used [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Gill | Institute of Mathematics and Its Applications[END_REF]. We will present here a preconditioned and re-orthogonalized implementation of this algorithm, used originally in [START_REF] Ben Dhia | Multimodeling of multi-alterated structures in the Arlequin framework[END_REF][START_REF] Elkhodja | Approches de structures complexes dans des cadres adaptés de la méthode Arlequin[END_REF]. The algorithm itself is presented in Alg. 1. From a numerical point of view, it consists on applying the CG method to the equation A + • Φ = b + , but modified in such a way that the constraint ( 15) is satisfied at each iteration. This can be done by choosing an initial solution guess, Φ (0) , that follows this constraint, and imposing that the descent directions p (k) are inside the null space of R I 2 . This guarantees that R I 2 Φ (k) = R I 2 Φ (0) , and thus Eq. ( 15) is satisfied. A suitable choice of Φ (0) is

Φ (0) = -R I 2 R I 2 T R I 2 -1 R 2 T • F 2 , (21)
and the descent direction condition can be achieved by applying an orthogonal projector operator,

Π R = I n RB 2 -R I 2 R I 2 T R I 2 -1 R I 2 T ( 22
)
on the vectors used to calculate the search directions.

After the constrained value of Φ is found, α can be calculated from the system's residual b + -A + Φ, using Eq. ( 16):

α = R I 2 T R I 2 -1 R I 2 T b + -A + Φ . ( 23
)
3.4 Implementation of the projected CG algorithm Some parts of the projected CG algorithm deserve a more detailed description. These include the numerical optimizations needed to avoid the explicit construction of the matrix operators A + and Π R , the effects of numerical fluctuations, a proper choice of a preconditioner, and the choice of the convergence criteria. We will focus now on these points. Indirect operator application: Two steps of the Alg. 1, involving the operators A + and Π R , at lines 7 and 11, are potential performance bottlenecks. These operators involve the evaluation of inverse (dense) matrices, and thus we want to avoid their explicit construction.

In the case of the system matrix A + , each application of a term C l K -1 l C T l can be exchanged by a multiplication by C T l , a call to the model l's solver, and a multiplication by C l -avoiding the explicit construction of the inverse operator. The same is true for the terms with the pseudo-inverse, but in this case one must use a version of the corresponding model's solver that takes the null space in consideration. For a linear elasticity model, the PETSc Krylov solvers can do so using the MatNullSpaceCreateRigidBody and the MatSetNullSpace functions, which create and attach the rigid body modes to a given matrix, respectively. Another approach in the same context would consist for the user to indicate a set of additional Dirichlet boundary conditions blocking the rigid body modes while keeping the same amount of work on the rest of the displacement field.

For the projection operation, a similar decomposition onto a series of operations can be used. In this case, the n RB

2 × n RB 2 matrix R I 2 T R I
2 is small enough to be inverted explicitly: in the worst case for a 3D linear elasticity model, it is a (6 × 6) matrix. Again due to this small size, the inverted matrix can be stored and multiplied locally, avoiding a communication bottleneck when Alg. 1 is used in parallel. A similar procedure can Algorithm 1: ProjectedPCG: preconditioned and re-orthogonalized Conjugated Gradient algorithm, modified to project the search direction into the null space of the matrix R I 2 . The inputs are the linear system's matrix and right-had side, A + and b + (given here by Eqs. ( 17) and ( 19)). The preconditioner operator is represented by the matrix M P C , and the projection operator Π R is given by Eq. ( 22).

Input:

System matrix:

A + = C 1 K -1 1 C T 1 + C 2 K + 2 C T 2 System right-hand side: b = (C 1 u 0 1 -C 2 u 0 2 ) Output: Converged solution Φ k f /* Initialization */ 1 Initial solution: Φ (0) = -R I 2 R I 2 T R I 2 -1 R 2 T • F 2 ; 2 r (0) = b -A + • Φ (0) ; 3 z (0) = Π R • M P C • Π R • r (0) ; 4 ρ (0) = r (0) , z (0) 2 ; 5 p (0) = z (0) ; 6 for k = 0, 1, 2 . . . until convergence do 7 q (k) = A + • p (k) ; 8 γ (k) = ρ (k) / p (k) , q (k) 2 ; 9 Φ (k+1) = Φ (k) + γ (k) p (k) ; 10 r (k+1) = r (k) -γ (k) q (k) ;
/* Preconditioning and projection */ k+1) , z (k+1) 2 ; /* reorthogonalization of the descent */

11 z (k+1) = Π R • M P C • Π R • r (k+1) ; 12 ρ (k+1) = r (
13 for i = 0, 1, 2 . . . k do

14 β i = z (k+1) , q (i) 2 / p (i) , q (i) 2 ; 15 p (k+1) = z (k+1) -k i=0 β i p (i) ; 16
CheckConvergence();

be followed for Eqs. ( 21) and [START_REF] Juntunen | Nitsche's mthod for general boundary conditions[END_REF], where this inverse matrix also appears.

Effects of numerical fluctuations:

A classical CG algorithm relies on the fact that the descent directions p (i) are orthogonal to each other, under the dot product associated to the matrix A + . In the case of the projected CG algorithm presented in Alg. 1, another condition is imposed on these vectors: that they are inside the null space of R I 2 . This guarantees that the condition ( 15) is followed, and thus that the solution ( 14) exists.

Due to numerical fluctuations, though, these properties are not followed exactly. In the case of the first condition, this results only into a slower algorithm convergence. More importantly, in the case of the second con-dition, these numerical errors result in incorrect rigid body mode corrections (the term R 2 α in Eq. ( 14)), posing a serious problem to the robustness of the projected CG algorithm.

One way to avoid the first effect is to re-orthogonalize the descent directions (line 12 of Alg. 1). This change increases the memory cost of the algorithm, since we must save all the vectors q (i) . But, since system associated to the coupling space is the smallest one of the coupled problem, this should not be an issue. In the case of the second condition, the application of the projection operator directly on the residual (line 11 of Alg. 1), together with the reorthogonalization, guarantee by construction that the new descent direction p (k+1) is inside R I 2 's null space. Convergence condition: Usually, the convergence condition followed by a CG algorithm takes into account the relative or absolute convergence of the residual k+1) , by calculating its norm under the preconditioner dot product

r (k) = b + -A + Φ (
ρ (k) = r (k) , M P C r (k) , (24)
which tends to zero as we approach the exact solution.

Due to the projection operations of the projected CG algorithm, though, this norm does not tend to zero.

A more suitable choice of norm is the projected and preconditioned dot product, used in line 11 of Alg. 1:

ρ (k) = Π R r (k) , •M P C • Π R r (k) . (25)
Furthermore, while we might reach a convergence of the solution Φ (k) , this does not guarantee that the rigid body modes correction, R 2 α, has converged. The function CheckConvergence(), at the line 16 of Alg. 1, must then be altered to check the convergence of this correction term. Here, we chose a relative convergence test, comparing two consecutive values of the the corrections:

abs(|R 2 α (k+1) | 2 -|R 2 α (k) | 2 ) < RB |R 2 α (k+1) | 2 . ( 26
)
The effects of this choice of convergence check will be discussed in the section 6.1, together with the effects of the reorthogonalization. Preconditioner: Following the ref. [START_REF] Ben Dhia | Multimodeling of multi-alterated structures in the Arlequin framework[END_REF][START_REF] Elkhodja | Approches de structures complexes dans des cadres adaptés de la méthode Arlequin[END_REF], we used a preconditioner M P C based on the coupling matrix C m . Let us justify, from a mechanical point of view, this choice. Using Eq. ( 7), the matrix A + can be rewritten in the following manner:

A + = C m P 1 K -1 1 P T 1 + C m P 2 K + 2 P T 2 C m . ( 27
)
Recall that the coupling matrix C m is analogous to a stiffness matrix with rigidity similar to those of the materials studied, and note that the terms P 1 K -1 1 P T 1 and P 2 K + 2 P T 2 correspond to interpolations of the inverses of stiffness matrices into the coupling region. The matrix A + can then, in the context of choosing a preconditioner, be roughly approximated by C m .

The most straightforward preconditioner choice is the full inverse of the coupling matrix,

M P C = (C m ) -1 ,
but its efficient usage depends on solving yet another equation system inside the main for loop of the projected CG algorithm. This system is defined in the mediator mesh M, though, and in most cases it should be relatively cheap to solve it, when compared to the full model's systems. If this is not the case, another option is the usage of a Jacobi-like preconditioner, with

M P C = [diag (C m )] -1
. This choice avoids solving a new equation system, but is a rougher approximation. We tested and compared both choices (together with simulations without the preconditioner), and the results are presented in section 6.3.

Other domain decomposition methods

As we said previously, other domain decomposition methods can be used to solve the Arlequin problem, and the choice depends on the application and the method's strengths. One of such is the LATIN method, which is applied to the Arlequin problem in ref. [28]. This method has the advantage of not being susceptible to the rigid body modes problem seen for the FETI method, since its "decoupled" step changes the rigidity matrices (Eq. ( 33) in the reference). Due to this, it is compatible with external solvers that do not project out the rigid body modes. Furthermore, it is directly compatible with the proper generalized decomposition method [28]. On the other hand, due to these same changes to the rigidity matrices, the LATIN / Arlequin solver is arguably more intrusive: either the external solver must accept the (possibly dense) modifications of the rigidity matrices, either the external system must be solved internally using an algorithm similar to Alg. 1, using a A • x operation from the external solver.

Coupling matrix assembly and parallelization

Let us discuss now the parallelization of the discretized formulation construction, Eq. ( 5). The matrices K 1 , K 2 and the vectors F 1 , F 2 can be constructed in parallel using FEM libraries such as libMesh [START_REF] Kirk | libMesh: A C++ Library for Parallel Adaptive Mesh Refinement/Coarsening Simulations[END_REF]. We will present in this section the parallelized construction of the coupling matrices C l . From Eq. ( 6), we have that these matrices involve the integration of form functions associated to different, incompatible meshes, M and T l .

Figure 2: (color online) 2D example of an intersection between two elements E m ∈ M (blue) and E l ∈ T l (red), corresponding to a non-zero coupling term C l i,j . The intersection I (gray) does not follow the geometry of either of these two meshes, and thus it must be triangulated to allow the calculation of the coupling term.

We must, then, define a common mesh, over which we can project these form functions, to calculate the integrals using a Gaussian quadrature. This common mesh, which we will note I m , must be constructed before the proper assembly of the coupling matrices.

The form functions ν m i and ν l j used in Eq. ( 6) are defined on the elements E m ∈ M and E l ∈ T l , respectively. For two meshes formed by three dimensional elements, the coupling terms C l i,j will be non-zero only if their intersection defines a polyhedral region I (i.e., it will be zero if the intersection is empty or if it defines a lower dimensionality entity, such as a point or a surface). The integral Eq. ( 6) must be calculated over this intersection, but in the general case it does not follow the element geometry of either M or T l (see Fig. 2 for a 2D example), and we cannot calculate the integrals using directly the original form functions.

We need, then, to decompose each region I into a form that can be treated using a Gaussian quadrature. This can be done by triangulating each intersection region I, defining an element set {E I }, Fig. 2. The union of all these intersection triangulations defines an intersection mesh, I m . Each element E I is associated to a set of quadrature points {q p } and a form function basis {ν I k }. The latter can be used to rewrite the form functions ν m i and ν l j . For linear form functions, this results Figure 3: 2D example of the projection of the evaluation of a form function ν m 1 of the element E m at a quadrature point q p of the intersection element E I . In the case represented, ν m 1 (q p ) is equal to the form function ν I 2 of the intersection elements, evaluated at q p and rescaled by a constant λ.

in linear combinations,

ν m i (q p ) = k λ i,m k ν I k (q p ), ν l j (q p ) = k λ j,l k ν I k (q p ), (28)
where the weights λ depend on the positions of the intersection elements' quadrature points inside the elements E m and E l . See Fig. 3 for a visualization of Eq. ( 28) for 2D elements. The coupling terms C l i,j can be rewritten by inserting the linear decompositions of ν m i and ν l j into Eq. ( 6). The result is a sum of integrals:

C l i,j = I k,k λ i,m k λ j,l k C I k,k , C I k,k = Ω c 12 κ ε ν I k : ε ν I k + 1 e 2 ν I k • ν I k dΩ. ( 29
)
Each C I k,k corresponds to an integral involving the form functions associated to an intersection element E I , and is defined on its region. It can, then, be calculated using its quadrature points. Since these integrals are associated to the intersection mesh I m , it is reasonable to use its partitioning over all the processors to distribute the calculations. Another possibility would be to use the partitioning of either M and T l , which determines the distribution of C l 's degrees of freedom, but doing so would leave idle processors, since not all of their parts are involved in the coupling. Alg. 2 shows a possible implementation of the parallelized coupling matrix assembly. It takes as input three meshes, the system mesh T l , the mediator mesh M and the intersection mesh I m , and an intersection table iT bl. This table encodes the relations between these three meshes, associating to each element E I ∈ I m the pair of parent intersecting elements {E m , E l }, E m ∈ M and E l ∈ T l . The parallelized loop on line 5 of Alg. 2 uses this information to calculate the contribution of E I using Eq. ( 29). Finally, the output is the parallel coupling matrix C l , with its row and column degrees of freedom following, respectively, the partitioning of M and T l . For simplicity sake, we suppose in Alg. 2 that all the processors have the data of all the elements E m and E l . Some FEM libraries do this by default to reduce the number of communications. If this becomes too memory intensive, reduced versions of M and T l , containing only elements that intersect the Ω c 12 coupling region, can be used.

Algorithm 2: BuildCoupling : Parallelized construction of the coupling matrix C l , following Eq. ( 29). The inputs are the system mesh T l , the mediator mesh M, the intersection mesh I m , and the intersection table iT bl. The latter contains, for each element E I ∈ I m the intersecting pair of elements {E m , E l }, with E m ∈ M and E l ∈ T l . The output is the coupling matrix C l .

Input:

System mesh T l , mediator mesh M Intersection mesh I m , intersection table iT bl

Output:

Coupling matrix C l 1 Partition the meshes T l , M;

2 Partition the mesh I m and the intersection table iT bl; The construction of each operator C l depends on building all the intersections between the two model's meshes.

In general, it is useful to separate this build process into two parts: a search step, which returns a list of intersecting elements, and a build step, which takes this list and build the geometrical intersection. The latter step is essentially the same, independently from the choice of search algorithm. Several serial search algorithms for two meshes exist in the literature [START_REF] Samet | Hierarchical spatial data structures[END_REF][26] [START_REF] Zomorodian | Fast software for box intersections[END_REF][18], but we must take into account the fact that our intersections are restricted to the coupling region Ω c 12 . This region can be represented by a third mesh, T C , and the problem can be formulated in a general form as: build all the intersections between three given meshes, T A , T B and T C . The mesh T C , though, does not have the same status as the other two, since it has no physical system associated to it. We will present in this section a parallel algorithm that takes this into account to solve the intersection search problem, but first we will describe some serial search methods, upon which our algorithm is based. Note that, as will be detailed below, the most critical component to guarantee the efficiency of the parallel algorithm is the distribution of the intersection construction workload over the processors.

Serial mesh intersection search

A straightforward method to find the intersections between two meshes would be to test each one of the element pairs (E A , E B ) for an intersection, resulting in an algorithm with complexity O(|T A | |T B |), where |•| is the number of elements of a mesh. More efficient algorithms organize the bounding boxes of the elements into hierarchical spatial data structures [START_REF] Samet | Hierarchical spatial data structures[END_REF] [START_REF] Massing | Efficient implementation of finite element methods on non-matching and overlapping meshes in 3d[END_REF], or use advancing front methods and neighbor information [START_REF] Gander | An Algorithm for Non-Matching Grid Projections with Linear Complexity[END_REF] to reduce the number of operations. In any case, the resulting set of n I polyhedrons must be checked for intersections with the third mesh T C . This step is effectively a second intersection search, and its complexity can be reduced by saving, during the first step, only the intersections found inside the region Ω c 12 . We will focus here mainly on the advancing front method from ref. [START_REF] Gander | An Algorithm for Non-Matching Grid Projections with Linear Complexity[END_REF]. This algorithm takes advantage of the neighbor structure of the meshes by noting that, if it is known that two elements A 1 ∈ T A and B 1 ∈ T B intersect, then A 1 's neighbors are likely to also intersect B 1 or B 1 's neighbors. An implementation adapted to include a restriction to a convex region Ω C (which can be, for example, the bounding box of T C ) is presented in Alg. 3. It works as follows: -Starting from an initial pair of elements intersecting inside Ω C , {A 1 , B 1 } (Fig. 4a), find all the intersections involving B 1 by recursively testing the elements around A 1 (Fig. 4b), taking care of marking which elements were tested already. -While doing this, test the new elements A i for intersections with B 1 's direct neighbors inside the region Ω C (line 15 of Alg. 3). Save the positive cases in a list of element pairs to be tested (Fig. 4c). This step and the previous one correspond to the inner loop of Alg. 3, between lines 10 and 15.

-After exhausting all the intersections possibilities involving B 1 , mark it as "already tested". -Repeat the algorithm with the next element pairs to be tested, ignoring the elements of T B already treated.

Notice that, in this algorithm, the meshes do not have the same standing: T B is used as a guide for the intersection search, while T A is only explored locally at each step. In most practical cases, the number of elements The implementation of the advancing front algorithm presented in Alg. 3 differs mainly on the dependency on the region Ω C . This is done to restrict the intersections to this region (line 13), and to stop the algorithm from exploring intersections outside of it (line 15). The second test might seem redundant at first, but it is needed because the advancing front algorithm is very sensitive on the direction over which the so-called front is advancing. Without it, is entirely possible that algorithm will add to the test list T estP airs a pair of elements (A 1 , B 1 ) that intersect outside the coupling region, and such that none of A 1 's neighbors intersect B 1 (see Fig. 5). In this case, B 1 will be marked as "already treated" after testing with A 2 , and the intersection with A 3 will not be found. This can be somewhat mitigated by expanding the test list at line 14 with knearest neighbors, but this slows down the algorithm, while not guaranteeing that this problem will always be avoided.

A i ∈ T A tested

Parallel mesh intersection search

Let us describe now the parallel intersection search algorithm, Alg. 4. As before, for simplicity we suppose that the geometrical information of all elements of T A and T B are accessible by all processors. The main idea behind this algorithm is to distribute the intersection search over the processors by dividing the meshes T A and T B into overlapping patches P A and P B . These patches are built in such a way that each overlapping patch pair (P A , P B ) covers exactly one element E C of the mesh T C (see Fig. 6). The patch construction process itself can be done using a simplified advancing front algorithm (Alg. 5), which is essentially the inner loop of 3, using the element E C as the guide and T i as the probed mesh. The intersection search is then reduced to an embarrassingly parallel problem, composed by |T C | completely independent problems of finding the intersections between (P A , P B ) inside E C . The only step in Alg. 4 which needs communications between the processors is the partitioning of T C -which is used to distribute the workload over the processors.

Algorithm 5: BuildPatch: patch construction algorithm based on ref. [START_REF] Gander | An Algorithm for Non-Matching Grid Projections with Linear Complexity[END_REF]. The inputs are the element E C which the patch will cover and the mesh T i from which it will be constructed. The list T estElem is a list of elements of T i to be tested for intersections with E C . The output is the patch P i .

Input:

Mesh T i ,

Scaling of the parallel mesh intersection algorithm

We will present now the implementation and scalability of our parallel intersection search and construction algorithm. In general, numerical geometry algorithms are susceptible to numerical imprecision and rounding errors, resulting into incorrect results, and this includes the intersection search and construction algorithms. The Computational Geometry Algorithms Library (CGAL) [30] avoids this problem by implementing its geometrical entities and algorithms using exact number types. We also chose it due to its Nef polyhedrons module. The Nef polyhedrons are defined by Boolean operations on a series of half-planes, and, by consequence, the intersection between two Nef polyhedrons can be constructed by applying an "AND" operation between them. This choice allows us to easily construct the intersection of an element triplet (E A , E B , E C ), regardless of the geometry of the elements chosen. The usage of exact geometry algorithms comes at the cost of more expensive operations, though, and as such a proper optimization and balancing of the code is needed. Exact geometrical tests, such as verifying if two elements intersect, are considerably cheaper than exact constructions of new geometrical entities, such as the construction of their intersection. Due to this, when using exact operations, an efficient intersection search algorithm is considerably cheaper numerically than the intersection construction itself. Table 1 illustrates this with the CPU time taken by our search algorithm and by the intersection construction, for meshes of the domains presented in Fig. 7. Due to this difference, we will study the scaling of the search algorithms and of the intersection construction separately. Also, we will present a mesh repartitioning based on the number of intersections that allows us to properly distribute the workload of the intersection construction. All simulations presented here were done on the FUSION cluster, using Intel Xeon E5-2670v3 @ 2.30 GHz.

Intersection search:

154 s Intersection build: 46351 s

Table 1: CPU time of the intersection search and construction steps, for the meshes described in Fig. 7 and on Table 2.

Mesh Elements Intersections

T 7, and for the resulting intersection mesh. For the latter, the number of intersections is also shown.

Intersection search scaling: Strong scaling is the property of an algorithm to solve the same problem n times faster when using n times as many processors. We tested the strong scaling of the intersection search algorithm Alg. 4 using meshes of the domains presented in Fig. 7. The two system meshes, T 1 and T 2 , have ∼ 1.75•10 6 and ∼ 4.38•10 5 elements, respectively, and the coupling mesh I C has 2.7 • 10 3 elements. The resulting intersection mesh (not pictured) has 6.66•10 6 elements, which is considerably larger than the other three meshes used to generate it.

We tested the algorithm using two different search methods for the intersections between the patches P A and P B : with the "naïve", direct pair search method (testing all the the element pairs for intersections), and with the advancing front method, described in Alg. 3.

Figs. 8a presents the wall time vs. number of CPUs for these methods (respectively, the red lines with the crosses, and the blue lines with circles) and the dashed black lines show the ideal, O(n) complexity scaling. The wall time represents, effectively, the time taken by the slowest processor, and the difference between it and the time taken by the fastest processor in each simulation (represented by the thin lines in these figures) gives to us an idea of how the workload is distributed in each simulation -and hence of the implementation's efficiency.

The small differences between the wall times and the fastest processor times when compared to the O(n) complexity curve indicate that the algorithm is scalable. This reflects the fact that the processors do not communicate during the parallelized loop in Alg. 4, as we discussed in sub-section 5.2. These curves diverge from the O(n) line, though, as the number of processors increase. This is due to an non-ideal distribution of the search workload, which is defined by I C 's partitioning. Here, we used libMesh's default partitioning weights, which associates to each element a weight equal to its number of nodes. This choice of weights has no infor-mation about the intersections associated to element E C ∈ I C , and this leads to an uneven workload distribution as the number of processors (and hence partitions) increases. This diverge is small enough, though, that it poses no problem given the speed of the intersection searches.

Notice that the advancing front algorithm is more than an order of magnitude slower than the direct pair search method. This happens because the advancing front algorithm needs to construct an intersection to stop it from exploring intersections outside of the coupling region (as was discussed in sub-section 5.1), which increases the algorithm's cost.

Intersection construction scaling:

The blue curves with circle markers in Fig. 8b present the wall time of the intersection construction as a function of the number of processors, for the same meshes as above, and using libMesh's default partitioning weights. This graph presents the same behavior as Figs. 8a, but since the intersection construction takes more time than the Table 3: Wall and fastest processor times for the intersection construction. The latter is least time taken by a processor during the simulation. Done with 336 processors.

search, the uneven workload distribution becomes more pronounced.

Differently from the intersection search, though, we know at this step how many intersections are associated to each element E C ∈ I C . The time taken to process each element E C is directly proportional to this quantity, and hence using it as the weights of the partitioning results in a more balanced construction algorithm. This is illustrated by the red curves with plus markers in Fig. 8b. While there is still some unbalancing for a large number of processors when using this new choice of weights, it is considerably smaller than before. This results in an algorithm that is ∼ 20% faster with the repartition (see Table 3, with the intersection constructions times for 336 processors.). More specifically, the wall time for the intersection construction with repartition is 2195 s on 24 processors and 152 s on 336 processors, while it is respectively 2182 s and 188 s without repartition. The speedup is therefore 14.4/14 > 1 with repartition and 11.6/14 ≈ 0.83 without repartition. Finally, it can be concluded that the strategy proposed for the intersection search and construction makes the coupling algorithm strongly scalable on 24 to 336 processors.

Numerical results

Let us now present our numerical results. We will study the properties and scaling of the FETI / Arlequin solver, based on the projected Conjugated Gradient (CG) algorithm (Alg. 1), and present an application case composed by a macroscopic, homogeneous model coupled with a microscopic heterogeneous model representing a polycrystalline material. To implement the Arlequin and FETI algorithms, we used the libMesh library [START_REF] Kirk | libMesh: A C++ Library for Parallel Adaptive Mesh Refinement/Coarsening Simulations[END_REF]. All simulations presented here were done on the FU-SION cluster, using Intel Xeon E5-2670v3 @ 2.30 GHz. The libraries that we developed to do these simulations (including the ones used for the intersection search scaling in the previous section) are freely available at the "Code Arlequin" GitHub repository [START_REF]Code Arlequin[END_REF].

Weak scalability is the property of an algorithm to require a constant wall time to solve a problem with n times as many degrees of freedom over n times as many processors. We will start by studying the weak scaling of this algorithm for different implementations -with and without the reorthogonalization operation, and for different convergence tests. After this, we will analyze these implementations in more details, by focusing on the behavior of a system where their differences are obvious. Finally, we will study the computational impact of the different preconditioner choices.

For all these cases, we applied the FETI solver on a coupled traction test, presented in Fig. 9. Both systems are modeled with the same 3D linear elasticity model, with a Young's modulus E = 200 GPa and a shear modulus µ = 80 GPa. Their domains Ω 1 and Ω 2 are discretized by the meshes T 1 and T 2 . Both models are coupled at the region Ω c 12 . A force density F = 100 kPa is applied on the face Γ F of the domain Ω 2 , on the x direction, while the face Γ 0 of the domain Ω 1 is clamped. Since the domain Ω 2 has no Dirichlet boundary conditions, we will use its rigid body modes for the projected CG algorithm. In all cases, the mesh T C of the coupling region has |T C | ∼ 3.5 • 10 3 elements.

Weak scaling

Usually, the weak scaling of iterative methods such as the CG algorithm is done by analyzing the evolution of the average time per iteration for different system sizes. This is done because increasing the system size increases the number of iterations of these solvers, and thus the wall time does not offer a suitable measurement of the scalability. In our case, though, , the duration of each iteration depends on the choice of solvers for the system models -which are assumed here to have weak scalings of their own. Thus, this quantity is not a suitable measurement for the weak scaling of solely the projected CG algorithm implemented here.

Considering this, we chose to use the number of iterations of the coupled solver until the convergence, n conv for the weak scaling. This quantity should stay stable if one of the model's meshes increases in number of elements while the mediator space stays the same. In the context of our traction test, this corresponds to increasing the number of elements of the mesh T 2 while keeping the mesh T 1 unchanged (and thus the mediator space). We will present in the following sub-sections the numerical costs of the coupling phase, compared to the model solver costs.

To do this scaling test, we used a mesh T 1 with |T 1 | ∼ 5 • 10 5 elements, and we varied the number of elements of T 2 from 7.2 • 10 3 to ∼ 4.9 • 10 5 , roughly doubling |T 2 | for each case. The number of processors was The simulations were done using the inverse coupling matrix preconditioner, M P C = (C m ) -1 . In all cases, we used the relative ρ (k) convergence check ρ(k) < ρ ρ(0), with ρ = 10 -8 , and, when applicable, the R 2 α convergence check from Eq. ( 26), with RB = 10 -6 . Fig. 10 presents n conv for each simulation. Fig. 10a shows the difference due to the reorthogonalization (or lack of thereof), while Fig. 10b shows the differences due to the different convergence check tests -namely, with only the ρ (k) check, only the R 2 α check, and with both checks. In Fig. 10a, both checks were used.

For the case without the reorthogonalization, the algorithm does not present any form of weak scaling. The number of iterations varies wildly between 365 and 866, with no correlation with the number of elements from the mesh T 2 . The reorthogonalization reduces and stabilizes considerably the number of iterations, which now vary between 52 and 66, which corresponds to a better weak scaling.

Concerning the convergence checks, using either only the ρ (k) check or the R 2 α check results in series of simulations with worse scaling than using both checks at the same time (notably for |T 2 | ∼ 1.7 • 10 4 ). While the norm of R 2 α is a stronger convergence parameter than ρ (k) (reflected by the larger number of iterations for the former), it does not supersede it, since the R 2 α check curve does not coincide with the curve with both checks. Considering this, we decided to keep both checks.

Reorthogonalization and convergence check

In this section, we will present in more details the results of the scaling for |T 2 | ∼ 1.7 • 10 4 . This series of simulations presents the largest difference of number of iterations until convergence between the different algorithm implementations (Fig. 10), and hence is better adapted to illustrate the differences between them. Fig. 11 shows the evolution of the projected residual ρ(k) (following Eq. ( 25)) and of the rigid body (RB) modes correction norm |R 2 α|, as a function of the it- eration. The black and the red curves show the results for reorthogonalized and non-reorthogonalized simulations, with both the convergence checks. The simulations with only one of the convergence checks follow the black line, ending after 8 iterations for the residual check, and after 35 for the RB check. Fig. 12 shows the final, deformed system for each implementation case, with deformations scaled by a factor of 10 3 .

We can see from Fig. 11 that the non-reorthogonalized implementation differs after only a few iterations from the other implementations, indicating that the orthogonality of the descents is lost shortly after the start of the simulation. After this point, we have recurrent peaks for the residual, with amplitude variations of a factor of ∼ 10 2 . The RB correction present smaller oscillations.

The reorthogonalization operation not only reduces considerably the number of iterations until convergence from 687 to 61 (for the same convergence checks), but it also eliminates these peaks from all but the first few iterations. As we can see from Fig. 12, the converged solution is similar in both cases. Using only one of the convergence checks yields results that are not properly converged. All these factors indicate that the reorthogonalization is needed to guarantee the robustness and efficiency of the projected CG algorithm, and that both convergence checks are needed to guarantee a correct result. to the preconditioner operations, and "Proj." to the projection operations, Π R .

Preconditioner effects

Before we pass to an applied case, let us discuss the effects of different preconditioners M P C on the projected CG algorithm convergence, and also analyze the (possibly) overhead time added to the algorithm for each choice. We tested two different cases:

M P C = C -1 m , based on Eq. (27); M P C = diag (C m ) -1
, which is a Jacobi-like approximation of the previous case. To use the first preconditioner efficiently, we must define a third equation system, using C m as the system matrix, and solve it. The Jacobi-like preconditioner, on the other hand, only adds a vector term-by-term multiplication.

To do these tests, we chose to use a coupled system similar to the previous traction test, with both system meshes T 1 and T 1 containing ∼ 1•10 6 elements. The simulations were done on 384 processors. The convergence parameters chosen were ρ = 10 -8 and RB = 10 -6 . Figs. [START_REF] Engquist | The Heterognous Multiscale Methods[END_REF] show the residual ρ (k) and the RB correction |R 2 α| 2 for all simulation cases. For reference, these figures also show the results with no preconditioning and ρ = 10 -8 and RB = 10 -4 (this case did not converge for RB = 10 -6 ). As one should expect, the preconditioned simulations present much lower values of ρ (k) , and the R 2 α correction norm goes down rapidly to reasonable values for the preconditioned simulations, while it stays slightly above the correct values during the non-preconditioned simulation.

This can be explained by the imprecisions of the solution approximation Φ (k) : Eq. ( 23) has been deduced for the real solution Φ, and it will only result on valid rigid body modes corrections if Φ (k) is a good approximation of it -which is false during the first iterations of the solver in all cases. Since the non-preconditioned simulation takes longer to yield a reasonable approximation of Φ, it presents the largest interval with incorrect R 2 α corrections -and hence the larger oscillations and slower convergence. Table 4 presents the timing and iteration statistics for each of the preconditioned simulations, including the timing of their main components: the system solvers, the preconditioner, and the projection operators. In both cases, the projection operation amount to less than 0.4% of the wall time. This indicates that our implementation of the projection operator Π R does not incur a performance bottleneck, and that the algorithm optimization should focus on reducing the cost of the models' solvers and the number of iterations of the coupled solver. We should stress here that the costs of model solvers and of the M P C = C -1 m preconditioner depend strongly on the number of processors chosen for the models and the coupling. Note that the matrix C m is inverted using the generalized minimal residual method (GMRES), PETSc's default parallel iterative linear equation system solver.

Regarding the preconditioner costs, the M P C = C -1 m case is considerably more expensive than the Jacobilike M P C = diag (C m )

-1 case (and even more expensive than the solves of the models). This is compensated by the reduced number of iterations, and hence reduce number of calls to both models' solvers, resulting into a simulation ∼ 19.3% faster in the former case. The high cost of the M P C = C -1 m preconditioner is mainly due to the relatively low sparsity of this matrix when compared to the other matrices, which reduces the efficiency of the linear solver associated to it. If such costs becomes too high, the results for the Jacobi-like preconditioner M P C = diag (C m )

-1 show that it is a viable numerical alternative.

Application: crack test

For an application test, we considered the coupled crack test, with the meshes shown in Fig. 14. The macroscopic system is described by an homogeneous 3D linear elasticity model, with a Young's modulus E = 200 GPa and a shear modulus µ = 80 GPa. The corresponding domain Ω 1 is represented by the mesh T 1 (Fig. 14a), with |T 1 | ∼ 3.4 • 10 4 elements. Its leftmost side is clamped, while displacements are imposed on right side holes, on the ẑ and -ẑ directions. The microscopic system is described by an heterogeneous 3D anisotropic linear elasticity model, with physical parameters c 11 = 198 GPa, c 12 = 125 GPa, c 44 = 122 GPa. Its domain Ω 2 (Fig. 14b) is divided into 250 crystals, each with its own random anisotropy direction following an uniform direction distribution on a sphere. The corresponding mesh has |T 2 | ∼ 7.2 • 10 6 elements. This model is located at the crack junction of the macroscopic mesh T 1 , and both meshes are coupled by the region marked in Fig. 14c, with a coupling mesh with |T C | ∼ 1.1 • 10 3 elements. Overall, the domains were meshed in such a way that the mesh T 1 's element size is ∼ 10× smaller near the crack junction than on the outer region, while the mesh T 2 's element size is ∼ 10× smaller than the smallest element size of T 1 . For the preconditioner, we used M P C = C -1 m . The simulation was run on 96 processors, and converged after 63 iterations, with the convergence parameters, we used ρ = 2 • 10 -4 and RB = 10 -4 . Fig. 15 present the resulting deformed meshes, as well as the von Mises stress for both the models. The displacement conditions on the macroscopic mesh T 1 result on a coherent deformation of the mesh T 2 , which presents a maximum von Mises stress near the tip of the crack. We can also visualize the grain anisotropy effects near this region, with different grains presenting different deformation responses.

Conclusion

We presented in this article a fully scalable framework to assemble and solve a coupled problem using the Arlequin method. The two main components of this framework are a new parallel algorithm to find the intersections between the meshes of the coupled models and a more detailed and scalable implementation of a FETI / Arlequin solver. While the total CPU time costs of the intersection construction algorithm are considerable (mainly due to the choice of the exact geometry libraries), this algorithm follows a strong scaling, which guarantees a fast assemble of the coupling operations if a suitable number of processors is chosen. Concerning the FETI / Arlequin solver, we presented a parallel implementation which has a weak scaling and for which most of time cost is due to the calls of the models' solvers. This translates into an algorithm with a small numerical overhead for the coupling operations, and which will not worsen the scalability of the models' solvers. This is possible due to our detailed study of the optimizations needed to avoid the numerical bottlenecks of the algorithm and a careful choice of convergence criteria. The usage of the FETI method also allows this algorithm to be easily adapted to multi-physics problems with very different coupled model pairs, such as linear / non-linear models, atomistic / continuum and deterministic / stochastic model pairs. Furthermore, the formulation of the coupling method is completely independent from the models, referencing them only through the calls to their solvers. This allows the usage of each models' own solver and code implementation without changing the coupling algorithm. The deterministic / stochastic model pair application is especially important for the numerical homogenization of random materials [START_REF] Cottereau | A stochastic-deterministic coupling method for continuum mechanics[END_REF][START_REF] Cottereau | Numerical strategy for unbiased homogenization of random materials[END_REF][START_REF] Cottereau | A Stochastic-deterministic Coupling Method for Multiscale Problems. Application to Numerical Homogenization of Random Materials[END_REF]. Another interesting extensions of this algorithm consist on generalizing it for a parallel solver with time coupled models, or for reduced-order solvers, such as the PGD [28]. Our future works will focus on these applications.

Finally, the work presented here can be extended to other physical couplings using the formulation from Eq. (1), such as the ones presented in refs. [START_REF] Wellmann | A two-scale model of granular materials[END_REF][START_REF] Chamoin | Ghost forces and spurious effects in atomic-to-continuum coupling methods by the Arlequin approach[END_REF][START_REF] Hu | Multi-scale nonlinear modelling of sandwich structures using the Arlequin method[END_REF], and the intersection search algorithm can be adapted to other numerical methods, such as the Nitsche method [START_REF] Massing | Efficient implementation of finite element methods on non-matching and overlapping meshes in 3d[END_REF].

Figure 1 :

 1 Figure 1: (color online) Domains used in the Arlequin framework. The model domains are Ω 1 and Ω 2 , the overlapping domain is Ω 12 , and the coupling domain (marked in gray) is Ω c 12 .

35 oneach processor p do 6 7 foreach 8 {E 9 foreach

 6789 Associate C l rows to M's DoF's; 4 Associate C l columns to T l 's DoF's; Local parts of I m : I p m ; element E I ∈ I p m do /* Get the intersecting element pair */ m , E l } = iT bl[E I ]; DoF pair (i, j) of {E m , E l } do 10 Transform the form functions ν m i and ν l j to E I 's form function basis; 11 Calculate the contribution of E I ; 12 Add it to C l i,j ; 13 return C l 5 Mesh intersection search

Figure 4 :

 4 Figure 4: Intersection search: (a) algorithm starts at the intersection between A 1 and B 1 . (b) It explores the region around the initial intersection and finds all the intersections near it (gray area). (c) During this process, B 1 's neighbors (here, B 2 and B 3 ) are tested for intersections with the elements of T A being explored, and any pairs found are added to a list for the next iteration -in this case, {B 2 , A 4 } and {B 3 , A 7 } were added.

Figure 5 :

 5 Figure 5: (color online) Behavior of the advancing front algorithm without the intersection region check done at line 15 of Alg. 3. In this case, the intersection between B 1 and A 3 is ignored.

Figure 6 :

 6 Figure 6: (color online) Graphical representation of the patch construction and intersection used in Alg. 4, and described in Alg. 5.

Algorithm 4 :

 4 ParallelInterSearch: Parallelized intersection search. The inputs are the meshes T A , T B and T C . The latter mesh is used to partition the search over the processors. The function InterSearch can be any serial search function, adapted to save only intersections inside a a given convex region (such as Alg. 3). The output is a table iP airs, containing the intersecting pairs.

4 6 7 BuildPatch 9 InterSearch

 4679 Input: Meshes T A , T B and T C Output: Intersection pairs table iP airs 1 Partition the mesh T C ; 2 oneach processor p do 3 Local partition of T C : T p C ; Local intersection pairs: iP airs p ; 5 foreach element E C ∈ T p C do Overlapping patches: P A , P B ; /* Build the patches */ (T A , E C , P A ); 8 BuildPatch (T B , E C , P B ); /* Search the intersections */ (P A , P B , E C , iP airs p ); /* Gather intersections lists */ 10 Gather(iP airs p , iP airs); 11 return iP airs

Figure 7 :

 7 Figure 7: (color online) Domains of the meshes used for the intersection search and construction scaling: (a) macroscopic and (b) microscopic domains, and (c) the coupling region. The black lines indicate their geometrical position compared to one another. The number of elements of the corresponding meshes are shown in Table2.

Figure 8 :

 8 Figure 8: (color online) Intersection algorithm strong scaling, wall time vs. number of processors. (a) Scaling of intersection search, for different methods used to find the intersections between two patches P A and P B : test all element pairs and advancing front method (b) Scaling of intersection construction, with and without the repartitioning step. The wall times are indicated by the thick lines with symbols. The thin, dashed lines with symbols, marked with "Min. time" in the legends, indicate the time taken by the fastest processor in each simulation.

Figure 9 :

 9 Figure 9: (color online) Domains of the meshes used for the traction coupling tests: (a) schematic representation of the test, (b) domain Ω 1 of the first model, (c) domain Ω 2 of the second model, and (d) the coupling region, Ω c 12 . The black lines on Figs (b), (c) and (d) indicate their geometrical position compared to one another.

Figure 10 :

 10 Figure 10: (color online) Weak scaling for different implementations of the projected CG algorithm: number of iterations until convergence vs. numbers of elements for the mesh T 2 : (a) with and without the reorthogonalization step, and (b) for different convergence tests ( only checking the residual convergence, only checking the R 2 α correction convergence (RB check), and doing both checks.). In (a), both checks were used. In all cases, |T 1 | ∼ 5•10 5 elements.

Figure 11 :

 11 Figure 11: (color online) Effect of the descent reorthogonalization on the evolution of the final (a) residual norm, ρ (k) , and (b) rigid body correction norm |R 2 α| 2 . The algorithm converged after 687 iterations without the reorthogonalization (red curves), and after 61 with it and both convergence checks (black curves). The main figures show the curves up to 200 iterations, while the insets show the full curves. The simulations with only the residual or the RB convergence check follow the black curves, stopping after 8 and 35 iterations, respectively. Number of elements in each mesh: |T 1 | ∼ 5 • 10 5 elements, |T 2 | ∼ 1.7 • 10 4 elements.

Figure 12 :

 12 Figure 12: (color online) Final, deformed traction test mesh for different projected CG algorithm implementations. The colors represent the von Mises stress. Number of elements in each mesh: |T 1 | ∼ 5 • 10 5 elements, |T 2 | ∼ 1.7 • 10 4 elements. The displacements were scaled by a factor of 10 3 .

Figure 13 :

 13 Figure 13: (color online) Effects of the preconditioner choice on (a) the residual norm, ρ (k) , and (b) the L 2 norm of the rigid body correction term, |R 2 α| 2 . Number of elements in each mesh: |T 1 | ∼ |T 2 | ∼ 1 • 10 6 elements. Preconditioned simulations were done with ρ = 10 -8 and RB = 10 -6 , and the non-preconditioned one with

Figure 14 :

 14 Figure 14: (color online) Meshes used for the crack test: (a) macroscopic homogeneous model mesh, T 1 , with the microscopic, polycrystalline and anisotropic model mesh inset, T 2 . (b) zoom on the microscopic mesh. (c) coupling region mesh, T C (black, thick lines). Number of elements in each mesh: |T 1 | ∼ 3.4 • 10 4 elements, |T 2 | ∼ 1.1 • 10 6 elements.

Figure 15 :

 15 Figure 15: (color online) Results of the crack test: (a) deformed macroscopic model mesh, T 1 , (b) deformed microscopic model mesh, T 2 (inserted into the macroscopic mesh), and (c) deformed microscopic model mesh with grains. In (a) and (b), the colors represent the von Mises stress (log scale) of the macro and micro models, respectively. Results obtained after 63 iterations, with ρ = 2 • 10 -4 and RB = 10 -4 , and using the coupling matrix preconditioner.

  ). (c) During this process, B 1 's neighbors (here, B 2 and B 3 ) are tested for intersections with the elements of T A being explored, and any pairs found are added to a list for the next iteration -in this case, {B 2 , A 4 } and {B 3 , A 7 } were added.Algorithm 3: AFIntersectionSearch: Intersection search algorithm based on ref.[START_REF] Gander | An Algorithm for Non-Matching Grid Projections with Linear Complexity[END_REF]. The inputs are the two meshes to be tested, T A and T B , and the region over which we want to restrict the intersections, Ω C . The list T estA is a list of elements of T A to be tested for intersections with the element B

		/* Test for intersection, and update
		iP airs if positive	*/
	13	UpdtInter(A t , B t , Ω C , iP airs);
		/* Add A t 's neighbors to test list	*/
	14	UpdtTest(A t , T estA, T reatedA) ;
		/* Test B t 's neighbors for new pairs */
	15	UpdtPairs(A t , B t , Ω C , N ewP airs, T reatedB);
		/* Update the pair test list	*/

t . The output is a table iP airs, containing the intersecting pairs {A t , B t }. Input: Meshes T A and T B , convex region Ω C Output: Intersection pairs table iP airs 1 List of pairs to test: T estP airs; 2 List of elements from T B already treated: T reatedB; 3 Search first intersection inside Ω C : {A 1 , B 1 }, A 1 ∈ T A , B 1 ∈ T B ; 4 T estP airs.insert({A 1 , B 1 }); 5 while T estP airs is not empty do 6 Elems. from T A to test: T estA; 7 List of new pairs: N ewP airs; 8 {A init , B t } = T estList.pop(); 9 T estA.insert(A init ); 10 while T estA is not empty do 11 Elements from T A already treated: T reatedA; 12 A t = T estA.pop(); 16 T estP airs.insert(N ewP airs); 17 T reatedB[B t ] = true; 18 return iP airs

  in Alg. 3 inner loop is considerably smaller than |T A |, resulting in a linear average complexity O(|T B |) [18]. The worst case scenario happens when all the elements of T A and T B intersect each other: in this extreme case we have |T A | |T B | intersections, and so O(|T A | |T B |) tests are done. From the memory point of view, no new data structures are created, and the algorithm uses O(|T A | + |T B |) space.

  Elements from T i to test: T estElem; 2 Elements from T i already treated: T reatedElem; 3 Search first intersection with E C : E init ∈ T i ; 4 P i .insert(E init );

		element E C
		Output:
		Patch P i
	1 5 T reatedElem.insert(E init );
		/* Add E init 's neighbors to test list	*/
	6 UpdtTest(E init , T estElem, T reatedElem);
	7 while T estElem is not empty do
	8	E t = T estElem.pop();
		/* Test for intersection, and update P i if
		positive	*/

9 UpdtPatch(E t , E C , P i ); /* Add E t 's neighbors to test list */ 10 UpdtTest(E t , T estElem, T reatedElem); 11 return P i

Table 2 :

 2 Number of elements of the meshes from Fig.

Table 2 .

 2

		10	2
											10	3
	Wall time (s)	10 10	1 0		Adv. front All pairs Min. time -Adv. front			Wall time (s)	10 2	Without repart. With repart. Min. time -Without repart.
						Min. time -All pairs						Min. time -With repart.
		10	-1	24	O(n)	48	96	192	384		24	O(n)	48	96	192	384
							Processors							Processors
							(a)							(b)

  2 norm of the rigid body correction term, |R 2 α| 2 . Number of elements in each mesh:|T 1 | ∼ |T 2 | ∼ 1 • 10 6 elements.Preconditioned simulations were done with ρ = 10 -8 and RB = 10 -6 , and the non-preconditioned one with ρ = 10 -8 and RB = 10 -4 .

	M P C	n conv Total time (s) Time / iter. (s) Model 1 (s) Model 2 (s) Precond. (s) Proj. (s)
	C -1 m	79	136.639	1.730	39.769	32.452	64.199	0.219
					29.105%	23.750%	46.985%	0.160%
	diag (C m ) -1	227	169.337	0.746	85.422	83.370	0.002	0.543
					50.445%	49.233%	0.001%	0.321%

Table 4 :

 4

Wall time costs for each component of the projected CG algorithm (Alg. 1), for different choices of preconditioner matrices, M P C . "Model 1" and model "Model 2" correspond to calls to the model's solvers, "Precond."

Table 5 :

 5 Wall times for the coupled crack test simulation, with convergence parameters ρ = 2 • 10 -4 and RB = 10 -4 . "System initialization" consists on the equation systems' initialization. The "Arlequin / FETI solver" time is subdivided into its main components, including the setup, the models' solves, preconditioning operations, and other (projections and vector manipulations)

	Operation	Time (s)
	System initializations	8.09
	Coupling matrix assemble 27.16
	I/O	52.48
	Arlequin / FETI solver	398.61
	Setup	2.53
	Solve model 1	52.99
	Solve model 2	340.14
	Precond.	2.05
	Other	0.90
	Total	486.34

Acknowledgements This work benefited from French state funding managed by the National Research Agency under project number ANR-14-CE07-0007 CouESt. The simulations were done using the CentraleSupelec-ENS Paris-Saclay's computing mesocenter, FUSION. We would like to thank its support team for the help provided during the simulations.


François Hild

Didier Elisabeth Amar

Stress Heterogeneity Effect on the Strength of Silicon Nitride

published or not. The documents may come   L'archive ouverte pluridisciplinaire