On the Efficient Reconstruction of Displacements in FETI Methods for Contact Problems

The final step in the solution of contact problems of elasticity by FETI-based domain decomposition methods is the reconstruction of displacements corresponding to the Lagrange multipliers for “gluing” of subdomains and non-penetration conditions. The rigid body component of the displacements is usually obtained by means of a well known but quite complex formula, the application of which requires reassembling and factorization of some large matrices. Here we propose a simple formula which is applicable to many variants of the FETI based algorithms for contact problems. The method takes a negligible time and avoids reassembling or factorization of any matrices.


Introduction
The FETI methods proposed by Farhat and Roux [1] turned out to be an efficient tool for the solution of large problems arising from the discretization of elliptic partial differential equations.Using FETI, the domain is partitioned into non-overlapping subdomains, an elliptic problem with Neumann boundary conditions is defined for each subdomain, and the inter-subdomain field continuity is enforced via Lagrange multipliers.The Lagrange multipliers are evaluated by solving a relatively well conditioned dual problem of a small size.Since the stiffness matrices of the subdomains are typically only positive semidefinite, the formulation of the dual problem enhances also some additional constraints which are associated with another set of variables that we shall call secondary multipliers.In the classical variants of FETI for linear problems or for small contact problems, the cost of the evaluation of the secondary multipliers is negligible.However, if the number of the subdomains that are used for the solution of contact problem is large, then it is not the case.
For the solution of problems the primal dimension of which is in billions, the evaluation of the secondary multipliers involves the products of matrices with some hundreds of millions columns.
The point of this note is to recall that the secondary multipliers can be obtained nearly for free if the constrained dual problem is solved by a mixed method such as the method of augmented Lagrangians.

TFETI for Contact Problems
Let us briefly describe a structure of the discretized contact problem without friction.Such a problem arises from the application of a variant of FETI under the assumption that the kernels of the subdomains are known a priori, as is always the case if the TFETI method [2] is applied.Let us denote the subdomain stiffness matrices K s , the matrices formed by their independent rigid body modes R s , and the subdomain load vectors f s , s = 1, . . ., N S .They can be assembled into global objects The non-penetration and the "gluing" of subdomains are described by the components B I and B E of the rectangular matrix The discretized primal problem with unknown displacements u then reads Using the standard notation where K † denotes a left generalized inverse of K so that KK † K = K, we can formulate the dual problem in Lagrange multipliers and Gλ = e. ( The problem of minimization on the subset of an affine space can be transformed into that on the subset of a vector space by means of the substitution λ = ν + λ, where λ satisfies G λ = e.The problem then reads Observing that the projector maps any ν which satisfies Gν = o to zero, we can use Q and P = I − Q to modify Eq. ( 6) to where ρ > 0 denotes a regularization parameter and T a nonsingular matrix which defines the orthonormalization of the rows of G.In what follows, we shall assume that GG T = LL T is the Cholesky decomposition of GG T and we can take T = L −1 .For convenience, we shall denote G = TG.The (GG T ) −1 action is called Coarse Problem (CP) solution.The action time and level of communication depend primarily on CP implementation, i.e. on the solution of GG T x = y.
If ν is a solution of Eq. ( 6) so that λ = ν + λ is a solution of Eq. ( 5), then corresponding displacements u can be evaluated by the formula where and the matrix B is formed by the rows of B that correspond to the active constraints, i.e., to the equality constraints and to those that satisfy λ I,i > 0.
If ≈ F , then the condition number of the Hessian of the cost function satisfies [3] where h and H denote the discretization and decomposition parameters, respectively.Let us recall Eq. ( 11) was the key ingredient of early proofs of the optimality of FETI (see Farhat, Mandel, and Roux [3] for linear problems and Dostal et al. [4] or references in [5] for contact problems).

Reconstruction Formula
Let the problem given by Eq. ( 8) be solved by SMALBE algorithm (see Dostal and Horak [6] or the book by Dostal [7]) which generates the approximations of the Lagrange multiplier µ for the equality constraints in the outer loop by means of the inexact solutions of auxiliary bound constrained quadratic programming problems in the inner loop by MPRGP (see Dostal and Schoeberl [8] or the book [7]).We get the solution ν of Eq. ( 8) and the corresponding multiplier µ for the equality constraints that satisfy the KKT condition Using Pν = ν, Qν = o, and simple manipulations, we get Let us denote so that After multiplication of Eq. ( 9) by B, we obtain and Let us rewrite Eq. ( 13), Eq. ( 15) and Eq. ( 17) If we sum the right hand sides of Eq. ( 17) and Eq. ( 18) and use P + Q = I, we obtain the right hand side of the Eq. ( 19).Thus Let us recall that G = L −1 G.Moreover, the nonzero entries of Bu correspond to λ I,i = 0, so the jump in the displacements can appear only on the contact interface, where λ I,i = 0.If we take into account components related to the sets E and I : λ I,i > 0, we get the final relation which is valid if Thus if we have λ and µ, we can evaluate α by Eq. ( 14) and α by Notice that the matrices on the right hand side of Eq. ( 14) are available throughout the solution procedure, while the evaluation of Eq. ( 10) requires effective assembling of G.If we work with the explicit G, i.e., with assembled L −1 G , then the final formula gets even simpler form The advantage of the latter formula is that its evaluation does not require the backward solve of the system with L, which can be useful when L is not supplied by the applied direct solution routine.

Numerical Experiments Using PERMON
As a model of 3D linear elasticity contact problem, we considered an elastic cube in contact with a rigid obstacle decomposed into 4096, 32768 and 110592 subdomains.
The numerical experiments have been computed by means of PERMON (Parallel, Efficient, Robust, Modular, Object-oriented, Numerical) [12], which is a collection of software libraries, uniquely combining QP (Per-monQP) and Domain Decomposition Methods (DDM) of FETI type (PermonFLLOP).Both modules are built on the top of PETSc [13], mainly its linear algebra part.
They were run on a Cray XC30 based supercomputer ARCHER [9] operated by EPCC.It consists of 4920 compute nodes.Each compute node contains two 2.7 GHz, 12-core Intel E5-2697 v2 (Ivy Bridge) processors and at least 64 GB of memory.Compute nodes are interconnected by the Aries interconnect using a Dragonfly topology.ARCHER's Rmax is 1642.5 TFlop/s in the Linpack benchmark.
There are several strategies for CP solution: iterative, direct solution, orthonormalization of G rows eliminating CP at all.Times required for standard computation of the amplitudes α of rigid body modes in the most efficient way using parallel direct solver SuperLU_DIST on the ARCHER supercomputer are depicted in Fig. 1.Although SuperLU_DIST was used for the CP solution and MPI subcommunicators were employed to decrease the communication cost of CP solution (see e.g.[10] and [11]), the CP solution becomes a bottleneck.In contrast to this, the simple substraction of two vectors takes a negligible time.

Conclusion
In this paper, we have introduced new efficient reconstruction formula for the computation of rigid body motions based on SMALBE multiplier avoiding the modified coarse problem matrix assembling, its factorization and solution.New formula results in significant time and memory savings which increase with the necessity to reconstruct the primal solution from the dual one repeatedly, e.g. for time dependent problems.

Fig. 1 :
Fig. 1: Times associated with CP solution determine costs for α computation on ARCHER.