On the efficient reconstruction of displacements in FETI methods for contact problems

MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE On the Efficient Reconstruction of Displacements in FETI Methods for Contact Problems David HORAK1,2, Zdenek DOSTAL1,2, Radim SOJKA2 1Department of Applied Mathematics, Faculty of Electrical Engineering and Computer Science, VSB–Technical University of Ostrava, 17. listopadu 15/2172, 708 00 Ostrava, Czech Republic 2IT4Innovations National Supercomputing Center, VSB–Technical University of Ostrava, 17. listopadu 15/2172, 708 00 Ostrava, Czech Republic da[email protected], [email protected], radim.sojk[email protected] DOI: 10.15598/aeee.v15i2.2322 Abstract. 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. Keywords Coarse problem, contact problems, displacement reconstruction, PERMON, SMALBE, TFETI. 1. 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. 2. 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 Ks, the matrices formed by their independent rigid body modes Rs, and the subdomain load vectors fs,s= 1, . . . , NS. They can be assembled into global objects K=diag(K1, . . . , KNS), R=diag(R1, . . . , RNS), f= [(f1)T, . . . , (fNS)T]T. (1) c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 237 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE The non-penetration and the “gluing” of subdomains are described by the components BIand BEof the rectangular matrix B= [B1,...,BNs],Bs=Bs I Bs E.(2) The discretized primal problem with unknown displacements uthen reads min 1 2uTKu −fTus.t. BIu≤o and BEu=o. (3) Using the standard notation F=BK†BT,G=RTBT, d=BK†f,e=RTf,(4) where K†denotes a left generalized inverse of Kso that KK†K=K, we can formulate the dual problem in Lagrange multipliers min 1 2λTFλ−λTds.t. λI≥o and Gλ= e. (5) 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 λ=ν+e λ, where e λsatisfies Ge λ=e. The problem then reads min 1 2νTFν−νT(d−Fe λ)s.t. νI≥ −e λI and Gν=o. (6) Observing that the projector Q=GT(GGT)−1G,(7) maps any νwhich satisfies Gν=oto zero, we can use Qand P=I−Qto modify Eq. (6) to min 1 2νT(PFP +ρQ)ν−νTP(d−Fe λ) s.t. νI≥ −e λI,TGν=o, (8) where ρ > 0denotes a regularization parameter and T a nonsingular matrix which defines the orthonormalization of the rows of G. In what follows, we shall assume that GGT=LLTis the Cholesky decomposition of GGTand we can take T=L−1. For convenience, we shall denote e G=TG. The (GGT)−1action is called Coarse Problem (CP) solution. The action time and level of communication depend primarily on CP implementation, i.e. on the solution of GGTx=y. If νis a solution of Eq. (6) so that λ=ν+e λis a solution of Eq. (5), then corresponding displacements ucan be evaluated by the formula u=K†(f−BTλ) + Rα,(9) where α=−(GGT)−1GBK†(f−BTλ), G=RTBT,(10) and the matrix Bis formed by the rows of Bthat correspond to the active constraints, i.e., to the equality constraints and to those that satisfy λI,i >0. If %≈ kFk, then the condition number of the Hessian of the cost function satisfies [3] κ(PFP +%Q|Im P)≤CH h,(11) where hand Hdenote 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). 3. 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 (PFP +ρQ)ν=P(d−Fe λ)−e GTµ.(12) Using Pν=ν,Qν=o, and simple manipulations, we get −e GTµ=P(F(Pν+e λ)−d) = P(Fλ−d).(13) Let us denote α=−(GGT)−1GBK†(f−BTλ),(14) so that GTα=−GT(GGT)−1GBK†(f−BTλ) = =Q(Fλ−d).(15) After multiplication of Eq. (9) by B, we obtain Bu =BK†(f−BTλ) + BRα= =d−Fλ+GTα,(16) and GTα−Bu =Fλ−d.(17) c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 238 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE Let us rewrite Eq. (13), Eq. (15) and Eq. (17) −e GTµ=P(Fλ−d),(18) GTα=Q(Fλ−d),(19) GTα−Bu =I(Fλ−d).(20) 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 −e GTµ+GTα=GTα−Bu.(21) Let us recall that e G=L−1G. 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 Eand I:λI,i >0, we get the final relation BRL−Tµ=BR(α−α),(22) which is valid if L−Tµ=α−α.(23) Thus if we have λand µ, we can evaluate αby Eq. (14) and αby α=α−L−Tµ.(24) 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 e G, i.e., with assembled L−1G, then the final formula gets even simpler form α=α−µ.(25) 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 Lis not supplied by the applied direct solution routine. 4. 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 (PermonQP) and Domain Decomposition Methods (DDM) 0.02 0.03 0.06 0.13 0.25 0.50 1.00 2.00 4.00 8.00 16.00 32.00 64.00 4096 110592 Time (s) (log scale) 32768 Number of subdomains (-) LU fact. of GGT(GGt)-1 action Assembling of GGT Fig. 1: Times associated with CP solution determine costs for αcomputation on ARCHER. 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 Grows 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. 5. 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. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 239 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE Acknowledgment The authors acknowledge the support of the Czech Science Foundation (GACR) project no. 15-18274S. The work was also supported by The Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPU II) project “IT4Innovations excellence in science - LQ1602” and from the Large Infrastructures for Research, Experimental Development and Innovations project “IT4Innovations National Supercomputing Center“ LM2015070”, and by the internal student grant competition project SP2017/169 “PERMON toolbox development III”. References [1] FARHAT, C. and F.-X. ROUX. A method of finite element tearing and interconnecting and its parallel solution algorithm. International Journal for Numerical Methods in Engineering. 1991, vol. 32, iss. 6, pp. 1205–1227. ISSN 1097-0207. DOI: 10.1002/nme.1620320604. [2] DOSTAL, Z., D. HORAK and R. KUCERA. Total FETI – an easier implementable variant of the FETI method for numerical solution of elliptic PDE. Communications in Numerical Methods in Engineering. 2006, vol. 22, iss. 12, pp. 1155–1162. ISSN 1069-8299. DOI: 10.1002/cnm.881. [3] FARHAT, C., J. MANDEL and F.-X. ROUX. Optimal convergence properties of the FETI domain decomposition method. Computer Methods in Applied Mechanics and Engineering. 1994, vol. 115, iss. 3–4, pp. 365–385. ISSN 0045-7825. DOI: 10.1016/0045-7825(94)90068-X. [4] DOSTAL, Z., T. KOZUBEK, V. VONDRAK, T. BRZOBOHATY and A. MARKOPOULOS. Scalable TFETI algorithm for the solution of multibody contact problems of elasticity. International Journal for Numerical Methods in Engineering. 2010, vol. 82, iss. 11, pp. 1384–1405. ISSN 00295981. DOI: 10.1002/nme.2807. [5] DOSTAL, Z., T. KOZUBEK, M. SADOWSKA and V. VONDRAK. Scalable Algorithms for Contact Problems. 1st ed. New York: Springer, 2017. ISBN 978-1-4939-6834-3. DOI: 10.1007/9781-4939-6834-3. [6] DOSTAL, Z. and D. HORAK. Theoretically Supported Scalable FETI for Numerical Solution of Variational Inequalities. SIAM Journal on Numerical Analysis. 2007, vol. 45, iss. 2, pp. 500–513. ISSN 0036-1429. DOI: 10.1137/050639454. [7] DOSTAL, Z. Optimal Quadratic Programming Algorithms, with Applications to Variational Inequalities. 1st ed. New York: Springer, 2009. ISBN 978-0-387-84806-8. DOI: 10.1007/b138610. [8] DOSTAL, Z. and J. SCHOEBERL. Minimizing Quadratic Functions Subject to Bound Constraints with the Rate of Convergence and Finite Termination. Computational Optimization and Applications. 2005, vol. 30, iss. 1, pp. 23–43. ISSN 0926-6003. DOI: 10.1023/B:COAP.0000049888.80264.25. [9] ARCHER. ARCHER [online]. 2017. Available at: http://www.archer.ac.uk/. [10] HAPLA, V., D. HORAK and M. MERTA. Use of Direct Solvers in TFETI Massively Parallel Implementation. In: 11th International Conference Applied Parallel and Scientific Computing (PARA 2012). Helsinki: Springer, 2012, pp. 192–205. ISBN 978-364236802-8. DOI: 10.1007/978-3-642-36803_14. [11] HAPLA, V. and D. HORAK. TFETI Coarse Space Projectors Parallelization Strategies. In: 9th International Conference Parallel Processing and Applied Mathematics (PPAM 2011). Torun: Springer, 2011, pp. 152–162. ISBN 978-3-64231463-6. DOI: 10.1007/978-3-642-31464-3_16. [12] PERMON team. Parallel, Efficient, Robust, Modular, Object-oriented, Numerical. PERMON [online]. 2017. Available at: http://permon.it4i.cz. [13] Portable, Extensible Toolkit for Scientific Computation. PETSc [online]. 2017. Available at: http://www.mcs.anl.gov/petsc/. About Authors David HORAK was born in Bohumin, Czech Republic. He has M.Sc. and Ph.D. degrees in Computer Science and Applied Mathematics and associate professor degree in Applied Mathematics at VSB–Technical University of Ostrava, 2007 Babuska price for Ph.D. thesis in Computer Science, since 2017 associate professor at Department of Applied Mathematics and since 2012 Researcher at IT4Innovations National Supercomputing Center, both at VSB– Technical University of Ostrava. His professional interest covers development and implementation of the scalable quadratic programming solvers using domain decomposition methods of FETI type, integral-discrete transforms. He is co-author of PERMON library. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 240 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE Zdenek DOSTAL was born in Olomouc, Czech Republic. He received his M.Sc. from Computer Science in 1970. He finished Ph.D. in 1979 in Mathematical Institute of the Czechoslovak Academy of Science, Prague, associate professor since 1992 from Charles University, Prague, full professor at VSB– Technical University of Ostrava since 1994, and Dr.Sc. from the Academy of Science of Czech Republic, Prague 2005. His research interests include numerical linear algebra, optimization, and computational mechanics. Radim SOJKA was born in Frydek-Mistek, Czech Republic. He received his M.Sc. from Comput. Mathematics in 2015 at VSB–Technical University of Ostrava. Since 2015, he is Research Asistant and Ph.D. student of Computer Science at IT4Innovations National Supercomputing Center, VSB–Technical University of Ostrava. His research interests include optimization, e.g. quadratic programming and integer programming for Vehicle Routing Problem and parallel implementation of optimization methods. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 241