É Alibart

S Picholle

G P Tanzilli

G L Puccioni

Lippi

T Wang

D Aktas

O Alibart

É Picholle

S Tanzilli

G P Puccioni

G L Lippi

Superthermal light from single-mode VCSEL

Superthermal light from single-mode VCSEL

Abstract: Light with superthermal properties is experimentally observed in the emission of a room-temperature, electrically-pumped single mode microcavity for pump values below continuous lasing. The observations are compared to Generalized Laser Rate Equations and stochastic simulations.

OCIS codes: 140.3945, 270.5290 Superthermal light has been recently recognized as a signature of Collective Effects (CE) in nanolasers since its origin is associated to the appearence of dipole-dipole correlations. The experimental observation of superthermal emission has been reported in the wings of emission pulses [START_REF] Jahnke | Giant photon bunching, superradiat pulse emission and excitation trapping in quantum-dot nanolasers[END_REF], accompanying the polarization switch [START_REF] Ch | Mode-switching induced super-thermal bunching in quantum-dot microlasers[END_REF][START_REF] Schlottmann | Exploring the Photon-Number Distribution of Bimodal Microlasers[END_REF], and in cw-pumped [START_REF] Kreinberg | Emission from quantum-dot high-beta micocavities: transition from spontaneous emission to lasing and the effects of superradiant emitter coupling[END_REF] microcavities. Thorough quantum-mechanical models, accompanied by numerical simulations, have also been proposed for systems with few emitters (and large coupling of the spontaneous emission into the lasing mode, i.e., β → 1) [START_REF] Jahnke | Giant photon bunching, superradiat pulse emission and excitation trapping in quantum-dot nanolasers[END_REF][START_REF] Kreinberg | Emission from quantum-dot high-beta micocavities: transition from spontaneous emission to lasing and the effects of superradiant emitter coupling[END_REF][START_REF] Leymann | Sub-and superradiance in nanolasers[END_REF]. A recent paper [START_REF] Protsenko | Collective Effects in Nanolasers Explained by Generalized Rate Equations[END_REF] proposes a set of Generalized Laser Rate Equations (GLREs) which includes the medium's polarization to characterize CEs in devices under more general conditions, i.e., independently of their β -factor, pump rate or relative relaxation rates (dipole vs. photon decay).

This approach is interesting since the standard adiabatic elimination, where the polarization variable is solved at equilibrium, is never entirely correct [START_REF] Oppo | Improved adiabatic elimination in laser equations[END_REF][START_REF] Oppo | Center-manifold reduction for laser equations with detuning[END_REF]. Thus, it stands to reason that even in those lasers -such as our microcavity -where the polarization decays much faster than the other variables (electromagnetic field and population inversion) a trace of its influence should survive. In this sense, the GLREs [START_REF] Protsenko | Collective Effects in Nanolasers Explained by Generalized Rate Equations[END_REF] represent a very interesting step forward for the description of small semiconductor-based devices (cf. [START_REF] Oppo | Center-manifold reduction for laser equations with detuning[END_REF] for a generalization of the simpler adiabatic elimination [START_REF] Oppo | Improved adiabatic elimination in laser equations[END_REF] to more realistic systems, including semiconductor lasers) and predict the occurrence of CEs in semiconductor lasers regardless of the limitations imposed by previous models [START_REF] Jahnke | Giant photon bunching, superradiat pulse emission and excitation trapping in quantum-dot nanolasers[END_REF][START_REF] Kreinberg | Emission from quantum-dot high-beta micocavities: transition from spontaneous emission to lasing and the effects of superradiant emitter coupling[END_REF][START_REF] Leymann | Sub-and superradiance in nanolasers[END_REF].

We report the experimental observation of superthermal light emitted by a single-mode microcavity (standard, electrically-pumped VCSEL) with a small fraction of spontaneous emission coupled into the lasing mode (β ≈ 10 -4 [START_REF] Wang | Dynamical Buildup of Lasing in Mesoscale Devices[END_REF]) in the transition region from spontaneous emission to a coherent, lasing field. The laser is operated at room temperature under cw operation with a very good control of all the parameters (pump, environment, etc.) and above threshold emits on a single polarization mode with rejection ratio up to 24dB (compared to the suppressed polarization). This strongly differentiates our observations from those already obtained in criogenically cooled, high-β microcavities [START_REF] Ch | Mode-switching induced super-thermal bunching in quantum-dot microlasers[END_REF][START_REF] Schlottmann | Exploring the Photon-Number Distribution of Bimodal Microlasers[END_REF], where superthermal light emission is accompanied by polarization hopping.

Fig. 1 shows a temporal trace of the emission collected from a VCSEL-980 (Thorlabs) VCSEL by a fast (9.5 GHz bandwidth, Thorlabs PDA8GS) detector, sampled by a digital oscilloscope (LeCroy Wavemaster 8600) at 0.1 ns rate (6 GHz analog bandwidth). It shows a sequence of (small, due to the detector sensitivity) independent pulsesbunched light -emitted without any particular regularity (caracterised by a nearly monotonically decreasing radiofrequency power spectrum -cf. Fig. 3, top left panel, of the Supplementary Information accompanying [START_REF] Wang | Dynamical Buildup of Lasing in Mesoscale Devices[END_REF]). The second-order, zero-delay autocorrelation g (2) (0) is obtained by equally splitting the laser output and detecting the coincidences with two single-photon detectors (APD idQuantique id100, ≈ 40 ps jitter) using and AND gate with a Time-to-Amplitude-Converter (Ortec 567, ≈ 15 ps timing resolution). The value obtained for the autocorrelation is g (2) (0) = 2.15 ± 0.05 [START_REF] Wang | Nontrivial photon statistics in small scale lasers[END_REF] in a regime where independent photon bunches are emitted, clearly indicating the superthermal nature of the emitted light.

Comparison to the theory [START_REF] Protsenko | Collective Effects in Nanolasers Explained by Generalized Rate Equations[END_REF] suggests compatibility with the influence of CEs in our system, since for the microcavity parameters (estimated to be K ≈ ×10 11 s -1 for the photon decay [START_REF] Wang | Dynamical Buildup of Lasing in Mesoscale Devices[END_REF], γ ⊥ ≈ ×10 13 s -1 for GaAs [START_REF] Coldren | Diode Lasers and Photonic Integrated Circuits[END_REF]) g (2) (0) values in excess of the thermal limit are expected (cf. Fig. 3 in [START_REF] Protsenko | Collective Effects in Nanolasers Explained by Generalized Rate Equations[END_REF]). The measured values are, however, larger than the limit estimated in [START_REF] Protsenko | Collective Effects in Nanolasers Explained by Generalized Rate Equations[END_REF] at g (2) (0) 2.01. On the other hand, such large values are compatible with stochastic simulations [START_REF] Puccioni | Stochastic Simulator for modeling the transition to lasing[END_REF] which, based on Einstein's semiclassical theory of stimulated emission [START_REF] Einstein | Zur Quantentheorie der Strahlung[END_REF] , do not include the dipole-dipole interactions. The stochastic simulations [START_REF] Puccioni | Stochastic Simulator for modeling the transition to lasing[END_REF] predict a dynamical output of the same kind as observed (Fig. 1), thus suggesting that the photon bunching be of dynamical nature. At the present time, however, it is possible that the CEs may also contribute to the dynamics [START_REF] Protsenko | Collective Effects in Nanolasers Explained by Generalized Rate Equations[END_REF], thus strengthening the superthermal emission properties. These points will be the object of additional investigations. 1 Currently at ... 2 Currently at ...

Fig. 1 .

 1 Fig. 1. Emission of bunched photons below the threshold for cw operation.


Paul Oumaziz

Pierre Gosselet

Pierre-Alain Boucard

Stéphane Guinard

A parallel non-invasive mixed domain decomposition -Implementation and applications to mechanical assemblies

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

It is a truism to say that industrialists are eager to conduct reliable simulations of their structures and therefore they need robust computational methods capable of handling models constituted by large amounts of degrees of freedom and involving complex nonlinear phenomena. Domain decomposition methods (DDM) are a natural way to distribute the computational load on several cores.

They also offer a framework where efficient iterative solvers can be implemented. In particular recent progress made it possible to design solvers with controlled convergence rate for linear symmetric positive definite systems for many DDM: see for instance [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF] for overlapping methods, [START_REF] Spillane | Automatic spectral coarse spaces for robust FETI and BDD algorithms[END_REF] for FETI and BDD methods, [START_REF] Klawonn | Adaptive Coarse Spaces for BDDC with a Transformation of Basis[END_REF][START_REF] Klawonn | Adaptive coarse spaces for FETI-DP in three dimensions[END_REF][START_REF] Beirão Da Veiga | Isogeometric BDDC preconditioners with deluxe scaling[END_REF] for FETI-DP and BDDC; also multipreconditioning proved to be an efficient alternative [START_REF] Bovet | Adaptive multipreconditioned FETI: scalability results and robustness assessment[END_REF] with possible application in more general contexts [START_REF] Bovet | Multipreconditioning for nonsymmetric problems: the case of orthomin and biCG[END_REF]. This paper focuses on problems involving frictionous contact between parts in small strain. Such a situation corresponds in particular to assemblies modeled at the scale of the connector (bolt....). For the domain decomposition simulation, a first possibility is to somehow decouple the contact from the domain decomposition by dealing with the contact inside subdomains, this was done in the alternate Schwarz case [START_REF] Badea | Schwarz method for dual contact problems[END_REF] or with the Balancing Domain Decomposition [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF]. In the view of dealing with assemblies, it seems more natural to aim at non-overlapping DDM able to deal with contact at the interfaces between subdomains. In this context, an important body of work is constituted by the contributions of Dostal and co-workers around the FETI(DP) method [START_REF] Dostál | Scalability and FETI based algorithm for large discretized variational inequalities[END_REF][START_REF] Dostál | Scalable FETI with optimal dual penalty for a variational inequality[END_REF][START_REF] Dostál | FETI based algorithms for contact problems: scalability, large displacements and 3D Coulomb friction[END_REF][START_REF] Dostál | Scalable TFETI with optional preconditioning by conjugate projector for transient frictionless contact problems of elasticity[END_REF], where the contact problem is rephrased as a constrained minimization problem dealt with a projected Krylov solver enriched by an efficient subdomain-based preconditioner and adapted coarse space. The alternative investigated in this paper is the Latin method [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF], which in the DDM language can be viewed as a non-overlapping optimized Schwarz method [START_REF] Lions | On the Schwarz alternating method. III: a variant for nonoverlapping subdomains[END_REF][START_REF] Gander | Optimized Schwarz Methods[END_REF], with dedicated treatment of the interface in order to insert complex mechanical behaviors.

Because they involve new practices and also because they require deep modification of the source code, DDM have not been massively adopted in commercial finite element software. In order to alleviate the latter problem, non-invasive implementations have been proposed, in particular since python driving of most software is possible. These implementations permit to demonstrate the methods on industrial problems and to exploit computational clusters, at the cost of reduced performance compared to true HPC implementations.

As an example, the global/local coupling [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Passieux | Local/global non-intrusive crack propagation simulation using a multigrid x-fem solver[END_REF][START_REF] Guguin | Nonintrusive coupling of 3D and 2D laminated composite models based on finite element 3D recovery[END_REF][START_REF] Guguin | On the computation of plate assemblies using realistic 3D joint model: a non-intrusive approach[END_REF] brings accuracy to the method of submodeling which is very popular amongst industrialists. It was applied as a complete solver in [START_REF] Duval | Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition[END_REF]. In fact the global/local coupling can be 1 viewed as a particular non-overlapping Schwarz method where the Robin condition is managed by a coarse representation of the domain [START_REF] Gosselet | Non-invasive global-local coupling as a Schwarz domain decomposition method: acceleration and generalization[END_REF]. This idea was exploited in [START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF] for the non-invasive implementation of the Latin method, the Robin conditions being simulated through the addition of elements on the interfaces of subdomains. In [START_REF] Oumaziz | A parallel noninvasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF], the extensibility of the non-invasive Latin method was studied by using the framework of [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF].

In this paper, we propose to demonstrate the efficiency of the non-invasive Latin method of [START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF][START_REF] Oumaziz | A parallel noninvasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF] with semi-industrial applications. Here semi-industrial refers to the fact that extremely fine models that would require HPC are out of the reach of our implementation. We restrict our study to infinitesimal strain and we make use of conforming meshes. Note that the Latin method can naturally be coupled with the mortar approach [START_REF] Bernardi | Domain Decomposition by the Mortar Element Method[END_REF][START_REF] Lacour | Two different approaches for matching nonconforming grids : the mortar element method and the FETI method[END_REF] in order to handle non-conforming interfaces, see [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF] for instance.

The paper is organized as follows: Section 2 presents the non-invasive Latin method; Section 3 explains the parallel implementation; finally Section 4 illustrates the implementation with several semi-industrial simulations involving many frictionous contact surfaces.

The substructured reference problem

We consider a family of non-overlapping subdomains (Ω E ) E∈[1,N ] of R 3 under quasi-static isotherm evolution and small perturbation assumptions. The subdomains are assumed to be made out of isotropic linear elastic material. Classical Lagrange conforming finite elements are used to approximate the subdomains' displacement field. For subdomain Ω E , we let K E be the stiffness matrix, u E be the vector of nodal displacement, f dE be the generalized forces associated to the given loads. We suppose that the degrees of freedom submitted to Dirichlet boundary conditions are treated by elimination.

Ω E Ω E Ω E f E E ,w E E f EE ,w EE f EE ,w EE f E E ,w E E f E E ,w E E f E E ,w E E Γ EE Γ E E Γ EE ∂Ω E E ∂Ω E E ∂ Ω E E ∂ Ω E E ∂ Ω E E ∂ Ω E E Figure 1:
Force and displacement fields defined on the edges of the subdomains

Subdomain Ω E interacts with its neighbor Ω E through the interface Γ EE = ∂Ω E ∩ ∂Ω E (Figure 1). The interface Γ EE is granted its own finite element discretization, and displacements and force distributions associated with its parent subdomains are defined: (w EE , w E E ) and (f EE , f E E ). These fields are connected to the subdomains by trace relations materialized by the operator N EE . We note E the set of all subdomains and G the set of all interfaces.

To describe complex load sequences and to handle the dependency of the solution with respect to history, a pseudotime is used. Thus nodal velocities ( ẇEE , uE ) are introduced with dot notations. These velocities are simply linked with displacement distributions through a backward Euler integration. If needed, the time step is indicated by a superscript. Typically:

ut+∆t E = u t+∆t E -u t E ∆t (1)
We use concatenated operators in order to easily manage many interfaces:

w E =     . . . w EE . . .     , f E =     . . . f EE . . .     , N E =     . . . N EE . . .     ,
where E spans all the neighbors of subdomain E

The mechanical problems to be solved on the subdomains can be written as:

∀Ω E ∈ E, ∀t, K E u t E = f d t E + N T E f t E w t E = N E u t E (3)
The mechanical behavior of the interfaces is given by a relation of the form:

∀Γ EE ∈ G,        f t EE + f t E E = 0, ∀t f t EE ∈ b EE ẇt E E -ẇt EE , t t; w 0 E E -w 0 EE , ∀t Initial conditions (4)
where the first relation is the balance of forces and the second relation represents the constitutive law of the mechanical behavior. The use of set notations for the behavior enables us to encompass most situations, among others perfect interfaces, contact, friction [START_REF] Blanze | A modular approach to structure assembly computations: application to contact problems[END_REF], cohesion [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF][START_REF] Saavedra | On a multiscale strategy and its optimization for the simulation of combined delamination and buckling[END_REF]. Block notations are used to simplify the writing of all the relations. As said earlier x E represents the gathering of all the (x EE ) E . x will represent the gathering of all the (x E ) defined on the subdomains, same procedure applies to operators:

x =     . . . x E . . .     , K =     . . . 0 K E 0 . . .    
, where E spans all subdomains.

(

We introduce operators which permit to make neighboring subdomains communicate: the operator A makes sums of interface vectors whereas the operator B makes differences:

(Af ) |Γ EE = f EE + f E E (Bw) |Γ EE = w E E -w EE (6)
These operators are orthogonal in the sense that ker(A) = range(B T ).

Thus, the discrete substructured problem can be written as:

∀t,      Ku t = f t d + N T f t Equilibrium of the subdomains w t = Nu t
Trace of the subdomain displacement

f t ∈ B T b(B ẇt , t t, Bw 0 ) Interfaces' behavior (7)
3 The non-invasive Latin method

Principle of the Latin method

So as to solve the substructured problem we use the Latin method [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF]. Two sets of interface distributions are defined :

A : (f t , ẇt ) t solutions to ∀t      Ku t = f t d + N T f t ẇt = N ut ut = (u t -u t-∆t )/∆t, + initial condition L : ( f t , ẇt ) t solutions to ∀t f t ∈ B T b(B ẇt , t t, B ẇt ) (8)
A groups the solutions to the linear equations set on the subdomains, A is an affine space called admissible space. L is a manifold containing the solutions verifying the behavior of the interfaces. In general the relations defining L are point-wise independent both in time and space.

From the initialization with s 0 ∈ A, an iterative scheme is used to reach the solution at the intersection of A and L. Partial solutions are alternatively searched in A and in L. This defines an iteration of the algorithm in two steps:

• the first step consists in searching an approximation s n = ẇn , f n in L starting from an approximation s n = (w n , f n ) in A. This stage is called local stage, it is computed by adding the relation of search direction:

f n -f n -k + V ẇn -ẇn = 0 (9)
• the second step is about finding an approximation s n+1 = (w n+1 , f n+1 ) in A starting from an approximation

s n = ẇn , f n in L.
This stage is called linear stage and it is computed by adding the relation of search direction:

f n+1 -f n + k - V ẇn+1 -ẇn = 0 (10)
k + V and k - V are two operators chosen by the user. k + V is chosen to be diagonal in order to benefit from the local properties of the relations defining L.

In order to ensure the convergence, a relaxation step is applied at the end of the linear stage:

s n+1 ← s n+1 + α(s n+1 -s n ) (11)
with 0 < α 1. The iterative method is represented graphically on Figure 2, where sn+1 is the result of the linear stage before relaxation. Proposition 1. The convergence of (s n+1 + s n )/2 is proved in [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF] when k - V = k + V are symmetric positive definite operators and α < 1 for maximal monotonic behaviors b.

A L k + V k - V s 0 s1 s 1 sn s n sn+1 s n+1 s th
Generally an indicator of error is used to quantify the energetic distance between A and L for two successive partial solutions. This indicator is defined as follows:

η = sn+1 -ŝ n 2 k - V sn+1 +ŝ n 2 k - V = ẇn+1 -ẇ T k - V ẇn+1 -ẇ + f n+1 -f T k - V -1 f n+1 -f ẇn+1 + ẇ T k - V ẇn+1 + ẇ + f n+1 + f T k - V -1 f n+1 + f (12)
3.2 Non-invasive multi-scale linear stage

Non-invasive linear stage

When combining the linear equations (8) of A together with the linear search direction [START_REF] Dostál | Scalability and FETI based algorithm for large discretized variational inequalities[END_REF], it comes out that solving the linear stage consists in solving independent Robin problems on the subdomains:

K + N T k - V ∆t N u t+∆t = f t+∆t d + N T f t+∆t + k - V ẇt+∆t + k - V ∆t w t (13)
Then compute

   ẇt+∆t = N ut+∆t = N(u t+∆t -u t )/∆t f t+∆t = f t+∆t + k - V ẇt+∆t -ẇt (14)
The search direction leads to a modification of the stiffness operator K by N T k -V ∆t N which is a non-standard term in industrial finite element software. In [START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF], we proposed a particular choice of the search direction to facilitate its implementation. The idea is to build the Robin extra stiffness by adding a layer of elements on the interface. This layer of elements will be also called a sole, it is illustrated on Figure 3. The intensity of the Robin parameter is adjusted by modifying the mechanical properties of the sole.

E 1 E 3 E 2 Γ 13 Γ 12

Multi-scale approach

In order to make the method scalable, long range effects need to be propagated, this is done via a multiscale extension [START_REF] Dureisseix | Une Approche Multi-échelles pour des Calculs de Structures sur Ordinateurs à Architecture Parallèle[END_REF][START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF][START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF]. We briefly recall the algorithm detailed in [START_REF] Oumaziz | A parallel noninvasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF]. The idea is to verify a weak form of the equilibrium of interface forces during the linear stage. Because of Saint-Venant's principle, a macro basis of displacements W containing rigid body motions is chosen, completed with simple deformation modes as extension or shearing modes. This new condition is written:

W T Af t+∆t = 0 ( 15
)
and it is taken into account in the search direction through the Lagrange's multiplier α:

f t+∆t -f t+∆t + k - V ẇt+∆t - ẇt+∆t + k - V A T Wα = 0 ( 16
)
Thus the problem to be solved at the linear stage consists in computing ẇt+∆t , f t+∆t knowing ( ẇt , f t ) and ẇt+∆t , f t+∆t verifying the following system:

Ku t+1 = f d + N T f t+∆t W T Af t+∆t = 0 w t+∆t = Nu t+1 f t+∆t -f t+∆t + k - V ẇt+∆t - ẇt+∆t + k - V A T Wα = 0 ẇt+∆t = w t+∆t -w t /∆t (17)
The algorithm to solve this system is fully detailed in [START_REF] Oumaziz | A parallel noninvasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF], beside extra communications and matrix operations during the initialization, it involves one extra global reduction per iteration.

Description of the local stage

Two types of interfaces are considered in our studies: perfect interfaces describing an interaction between subdomains of a same part, and contact interfaces describing physical interaction of subdomains of different parts. These two kinds of interfaces lead to two specific solving.

Perfect interface

For the perfect interfaces the problem consists in finding ẇt+∆t , f t+∆t knowing ẇt , f t and ẇt+∆t , f t+∆t verifying the system:

A f t+∆t = 0 Equilibrium of forces B w t+∆t = 0 Continuity of displacement f t+∆t -f t+∆t -k + V ẇt+∆t -ẇt+∆t = 0 Search direction ( 18
)
The solution is explicit and equal to:

ẇt+∆t = A T Ak + V A T -1 Ak + V ẇt+∆t -Af t+∆t f t+∆t = B T Bk + V -1 B T -1 Bk + V -1 f t+∆t -B ẇt+∆t ( 19
)
The computation of the velocity consists in solving a problem on two assembled soles with a Neumann conditions on the interface between them.

Contact interface

For the contact interface the problem consists in finding ẇt+∆t , f t+∆t knowing ẇt , f t and ẇt+∆t , f t+∆t verifying the system:

B n w t+∆t n + j n 0 Non interpenetration A f t+∆t = 0 Equilibirum of force B n f t+∆t n 0 Positive reaction force B τ f t+∆t τ = 0 If frictionless contact B n f t+∆t n T B n w t+∆t n + j n = 0 Signorini conditions f t+∆t -f t+∆t -k + V ẇt+∆t -ẇt+∆t = 0 Search direction w t+∆t = w t + ∆t ẇt+∆t , w 0 given Time integration ( 20
)
where j n is an initial gap. Choosing j n to be negative is one of the possibilities to simulate initial tension. We assume that the displacements and forces are written in a normal/tangential basis at each node of the interface, the n and τ subscripts refer to the normal and tangential components, respectively.

In the case of frictional contact the tangential condition B τ f t+∆t τ = 0 is replaced by the Coulomb relations which separate the sticking (first line) and sliding (second line) cases:

   f t+∆t τ < µ| f t+∆t n | then B τ ẇt+∆t τ = 0 f t+∆t τ = µ| f t+∆t n | then ∃λ > 0, B τ ẇt+∆t τ = -λ B τ B T τ -1 B τ f t+∆t τ (21)
where the use of normal (non bold) type refers to node-wise quantities, and µ is Coulomb's coefficient of friction. It is possible to define two node-wise indicators C n and G τ which permit to determine the status of nodes at the interfaces: C n gives information about the normal status (contact); in case of contact, G τ gives information about the tangential status (sticking/sliding). They are computed with known quantities and once they are determined the solution is explicit. The details of the solutions can be found in [START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF].

Parallel implementation 4.1 General environment of the developed tools

The chosen paradigm is to develop only with tools used by our industrial partners EDF and Airbus. In particular, we use code aster finite element software [START_REF]Code-aster[END_REF] and salome platform to generate and manage meshes within the MED format.

These two piecses of software are driven by python scripts permitting to develop a generic process, for the management of meshes, the mechanical computation and the post-processing.

In this part we describes the ideas of code-aster until the version 14. The further version may change the implementation of code-aster. The scripts are written in python with specific syntax to call the commands of code-aster. As the source of code-aster is written in fortran mostly, it includes a so-called supervisor (Figure 4) whose goal is to make the communications between the python and the fortran. The memory is shared with python. However, contrary to python, the aster objects are global in memory which lead to specific implementations that can be found in [START_REF] Oumaziz | Une méthode de décomposition de domaine mixte non-intrusive pour le calcul parallèle d'assemblages[END_REF].

Ideas of the developed tools

The tool we developed is composed of two parts:

• Preparation: starting from a monolithic mesh given by our industrial partners, the script detects all the group of nodes and elements needed for the further computations and distribute the data into individual files for each subdomain (data of the subdomain and its interfaces). If needed parts are subdivided into smaller subdomains. Remark. The recognition of the neighboring subdomains and the automatic construction of the interfaces and soles are crucial to have an automatic pre-processing phase. We have made the choice to detect the potential neighbors through the common boundary nodes of the subdomains. This easily permits to construct the interfaces and the associated soles. The automatic detection of neighbors for non-conforming interfaces is definitely not possible with that method and would require more evolved strategies.

An XML input file permits to define options and parameters both for the method and the structure. It is used to drive these two parts of the tools. This implementation is summed up in Figure 5.

Global mesh (MED format

Parallel processing

The non-invasive tool has been developed to launch parallel instances of python, each one driving a sequential version of the finite element software code aster [START_REF]Code-aster[END_REF]. The parallel instances of python are synchronized by the MPI protocol through the python package mpi4py [START_REF] Dalcín | MPI for Python: Performance improvements and MPI-2 extensions[END_REF]. Each MPI instance stands for a subdomain with its own interfaces. This principle is illustrated in Figure 6 with a simple four subdomain decomposition.

At the beginning of the computation, exchanges are required in order to assemble the coarse problem on each subdomain. Then each iteration of the Latin method is highly parallel and involves only one all-neighbors and one all-reduce communications:

1. The local is decomposed in two steps:

(a) An all-neighbor communication to exchange interface fields.

(b) A parallel computation to solve the local problems corresponding to its interfaces. Note that our implementation does not fully exploit the point-wise nature of the local equations and a finer degree of parallelism could be achieved (using one thread per interface node for instance).

2. The linear stage is decomposed in four steps:

(a) A parallel computation on each subdomain with the Robin condition.

(b) A global reduction to determine the macro lack-of-balance.

(c) A parallel computation to solve a coarse problem and to get the macro correction to be added to the displacement field.

Remark. In our case only sequential version of code-aster are used for each subdomain. One could think about calling parallel instance of code-aster for each subdomain, multi-threading or parallel solver included in code-aster would be a next step to improve the performance.

Automatic choice of search direction

The search direction, or in other word the Robin condition, is a parameter which strongly influences the convergence of the algorithm. It should represent the stiffness of the complement of the subdomain in the whole structure.

We propose a simple and automatic procedure to set this parameter for each interface. For parallel efficiency reason, our heuristic is only based on local data readily available to the processor and it does not require communication. Our method aims at having the interfaces develop a stiffness comparable to the subdomain's in their normal direction. Proposing more elaborate strategy will be the subject of future work.

Let us thus consider one interface of one subdomain. In our non-invasive implementation, the parameter of search direction is constructed by weighting the Young modulus of the sole whereas the Poisson coefficient is chosen to be the same as for the subdomain. During the preparation step, bounding boxes of soles and subdomains are determined (see Figure 7). Let ∆SD = [∆X, ∆Y, ∆Z] denote the bounding box of the subdomain, and ∆s = [∆x, ∆y, ∆z] the sole's. The mean normal to the interface is determined and denoted by n. The Young modulus of the sole E s is chosen with respect to the Young modulus of the subdomain E SD as follows:

E s |∆s • n| = E SD |∆SD • n| (22)
Figure 7: Example of bounding box with the subdomain in green and the sole in orange and the mean normal of the interface Remark. This automatic choice of weighting of the Young modulus does not work with all configurations and sometimes a manual choice is required. For example if the interface is circular, the mean normal is null. Moreover for warped subdomains, the bounding box may not represent the subdomain and the stiffness may not be approximated correctly.

Numerical examples

In this section, two numerical examples from industrial structures are presented. The first one corresponds to a simplified joint of a spatial launcher under two cycles of load in order to illustrate the accumulated sliding at the frictionous interfaces.

The second one is a joint which assembles the wings to the fuselage of an Airbus aircraft. With this example we illustrate the issue of overloaded subdomains which slow down the computation and we propose a simple load balancing strategy by re-decomposition.

Simplified joint of spatial launcher

The objective of this example is to illustrate the behavior of a joint composed of three bolts under two cycles of loadingunloading. It involves 30 subdomains linked through 53 interfaces (Figure 8). The three screws are decomposed in two parts linked with preloaded interfaces. This permits to ensure an initial state of compression between flanges. Contact interfaces are considered in the central part of the structure whereas perfect interfaces are chosen in the rest of the structure.

The structure is clamped on one side and a traction load is applied on the opposite side. The load follows two cycles as presented in Figure 9. 17 time steps are considered. The structure is discretized into 3 million degrees of freedom. All the parameters are given in Table 1.

The evolution of the error indicator is shown in Figure 10. Even though the multiscale approach speeds up the convergence at the first iteration, its effect decreases very quickly. The many non-linearities of frictional contact as well as the many timesteps slow down the convergence: more than 120 iterations are required to reach a criterion of 10 In Figure 11, we observe the preload of bolts at the first timestep. The screws are solicited with tensile stress whereas compression is observed in the surrounding plates.

Von Mises stress at the timesteps t 4 and t 12 (when the structure is the most loaded) is represented in Figure 12. Sliding 215 happens between plates and it causes shear stress in the screws, resulting in the structure being more solicited at the 12 th timestep than at the 4 th . Figure 13 presents the evolution of the longitudinal displacement at the end of the structure as a function of the load, we observe the accumulated sliding. During the loading phase, sliding appears and remains during the unloading phase.

Iterations Error indicator

Initial substructuring

This case permits to illustrate the feasibility of a resubstructuring in order to balance the size of subdomains. The structure was provided by Airbus. It represents the joint between the wings and the fuselage of an aircraft. The whole structure is composed of 47 subdomains and 89 interfaces, mostly frictionous. 10 bolts are considered. The structure is presented in Figure 14. It is clamped at the bottom and at the right. A traction load is imposed on the left. The screws and nuts are in titanium, whereas the rest of the structure is in aluminum. The load consists in one loading-unloading cycle discretized in 9 timesteps. All parameters are given in Table 2.

70 iterations are required to reach an error indicator around 10 -6 . We present in Figure 15 the stress of Von Mises for some timesteps. Once more, an accumulated sliding can be observed.

Re-substructuring of subdomains

Table 3 shows the number of nodes and the number of connected interfaces of the biggest subdomains. The biggest one represents almost a third of the whole structure and a quarter of the interfaces are connected to it. Thus the computation time is directly driven by this subdomain. For this computation, the total time is about 121.800 seconds. Therefore the objective is to re-decompose the large subdomains to reduce the number of degrees of freedom in each subdomain and also the number of interfaces per subdomain. The idea is to have new subdomains with as many dof as the smaller subdomains or to decompose the largest subdomains proportionally to their size in order to avoid the creation of too many subdomains. Thus the re-substructuring is limited by a user parameter M : the biggest subdomain is decomposed in at most M new ones, other subdomains are decomposed in proportion. This decomposition is realized through an existing function in code aster based on SCOTCH [START_REF] Pellegrini | Scotch: A software package for static mapping by dual recursive bipartitioning of process and architecture graphs[END_REF] whose only parameter is the number of subdomains.

In our case, we choose M = 10. The new structure is composed of 69 subdomains (Figure 16). With this new decomposition, we also ran 70 iterations and compared the CPU time. It permitted to reduce the CPU time to 18.450 seconds contrary to 121.800 seconds before, representing a gain of 6.6 from the initial substructuring. Figure 17 comparison of the error indicator between the two substructuring. The convergence of the first iterations are relatively similar. The multiscale approach grants scalability to the method. However the effect of the multiscale extension is decreased by the numerous non-linearities, and differences appear in term of convergence. to reach a plateau after 20 iterations. This is simply explained by the particular substructuring provided automatically by SCOTCH. For example in Figure 18, we show the new decomposition of the subdomain S 0 into 10 new substructures. One can remark that some subdomains are non-connected (composed of distinct parts).

The automatic weighting of the Young modulus cannot be efficient on such subdomains as the bounding box is not representative of the subdomain. Therefore, we decompose these problematic subdomains into new subdomains. With this, we obtain a decomposition of the whole structure into 72 subdomains. Such a decomposition is much more adapted to the automatic weighting of the search direction and we expect a better convergence. Therefore, the error indicator no longer reaches a plateau and the convergence follows the initial one but 6.6 times faster (Figure 17) In this paper we presented a multiscale non-invasive implementation of a mixed domain decomposition method. This implementation has been confronted to industrial structures with quasi-static and small strain assumptions, with many frictional contact interfaces. The non-invasive implementation of the mixed conditions is realized by adding of a layer of elements at the interface. We briefly reminded the non-invasive multiscale approach detailed in [START_REF] Oumaziz | A parallel noninvasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF] and the treatment of frictional contact [START_REF] Blanze | A modular approach to structure assembly computations: application to contact problems[END_REF][START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF].

260

Moreover we presented an automatic setting of the Robin parameter which is suited for non warped and homogeneous subdomains. Although this permits to avoid a manual parametrization of all search directions, we exhibited some limita-tions for certain shapes of subdomain. We also confronted the method to semi-industrial cases involving bolted assemblies with many frictional contact interfaces.

The next step will be to investigate the weighting of search direction in the case of heterogeneous materials in subdomains. The choice of only one Young modulus may not be efficient enough to ensure robust convergence. This needs also to be coupled with spectral coarse spaces such as proposed in GenEO approach [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF].

Figure 2 :

 2 Figure 2: Latin method -(s is the approximation before relaxation)

F

Figure 3 :

 3 Figure 3: Example of a subdomain with its soles

Figure 4 :

 4 Figure 4: Principle of the supervisor

Figure 5 :

 5 Figure 5: Description of the global implementation

Figure 6 :

 6 Figure 6: Principle of the implementation

Figure 10 :

 10 Figure 10: Evolution of the error indicator

Figure 11 :Figure 12 :

 1112 Figure 11: Preload of bolts -stress σ zz (MPa)

Figure 13 :

 13 Figure 13: Mean longitudinal displacement at the end of the structure

Figure 14 :

 14 Figure 14: Description of subdomains

Figure 15 :

 15 Figure 15: Stress of Von Mises (Pa) and deformed shape (×20)

Figure 16 :

 16 Figure 16: Re-decomposition of the structure

Figure 17 :Figure 18 :

 1718 Figure 17: Comparison of the error indicator

Table 1 :

 1 -5 . 0 t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 8 t 9 t 10 t 11 t 12 t 13 t 14 t 15 t 16 Parameters

	Figure 8: Description of the joint
	F d
	F max d
		t
	Figure 9: Load history
	Parameters	Values
	Young modulus E	100 GPa
	Poisson ratio ν	0.3
	F max d	10 MPa
	preload (left -middle -right) 0.3 -0.2 -0.4 mm
	∆t	1s
	dofs	3 millions

t

Table 2 :

 2 Parameters

	shows the

(a) Timestep t 2 (b) Timestep t 4 (c) Timestep t 6 (d) Timestep t 8

Table 3 :

 3 The new substructuring seems Size of the biggest subdomains

	Biggest subdomains Number of nodes Connected interfaces
	1 st	29,107	23
	2 nd	17,789	7
	3 rd	14,530	6
	4 th	11,030	10
	5 th	10,322	12
	6 th	6,861	6
	7 th	2,072	4

Núcleo Científico Multidisciplinario-DI, Universidad de Talca, Camino Los Niches km1, Curicó, Chile

LMT / ENS Cachan / CNRS / Université Paris-Saclay, 61 avenue du président Wilson, 94235 Cachan Cedex, France

Airbus SAS, CRT/XRXF,18 rue Marius Terce, 31300 Toulouse, FRANCE

Acknowledgments

We thank Airbus and EDF for their financial supports and technical help with code-aster. We also thank the project ECOS-CONICYT C17E04 for their financial support.