Extension of the parallel Sparse Matrix Vector Product (SpMV) for the implicit coupling of PDEs on non-matching meshes

Extension of the parallel Sparse Matrix Vector Product (SpMV) for the implicit coupling of PDEs on non-matching meshes G. Houzeauxa,∗, R. Borrella, J.C Cajasa,M.V´azqueza aBarcelona Supercomputing Center, Nexus II Building c/Jordi Girona, 29 08034 Barcelona (Spain) Abstract The Sparse Matrix Vector Product (SpMV) is one of the main operations of iterative solvers, and, in a parallel context, it is also the siege of point-topoint communications between the neighboring MPI processes. The parallel SpMV is built in such a way that it gives, up to round off errors, the same result as its sequential counterpart. In this regards, nodes on the interfaces (or halo nodes if halo are considered) are duplicated nodes of the same original mesh. It is therefore limited to matching meshes. In this work, we generalize the parallel SpMV to glue the solution of non-matching (non-conforming) meshes through the introduction of transmission matrices. This extension of the SpMV thus enables the implicit and parallel solution of partial differential equations on non-matching meshes, as well as the implicit coupling of multiphysics problems, such as fluid-structure interactions. The proposed method is developed similarly to classical parallelization techniques and can therefore ∗Corresponding author Email addresses: [email protected] (G. Houzeaux), [email protected] (R. Borrell), [email protected] (J.C Cajas), [email protected] (M. V´azquez) Preprint submitted to Computers & Fluids January 30, 2018 be implemented by modifying few subroutines of an already MPI-based code. According to the proposed framework, the classical parallelization technique appears as a particular case of this general setup. Keywords: Parallel sparse-matrix vector product, SpMV, MPI, Parallelization, Non-matching meshes, Coupling 1. Introduction At the algebraic level, historical methods to glue non-matching meshes are the mortar method [1, 2] and the Finite Element Tearing and Interconnecting (FETI) method [3], where the continuity of the solution is imposed through a Lagrange multiplier. See [4] for a comparison of both methods. The main drawback of these strategies is, apart from introducing additional unknowns to the original problem, their non-trivial implementations. The alternative method we propose was introduced in [5] and this work concentrates on implementation aspects, in a parallel computing environment. It is based on extending the parallel matrix-vector product through the introduction of transmission matrices, to express Dirichlet and Neumann couplings between non-matching subdomains. Also, these transmission can be built in such a way to obtain local and global conservation properties. The resulting method is implicit and can be implemented on the top of already existing parallelization methods for the sparse matrix-vector product (SpMV). Figure 1 shows some applications of the method, from the coupling of one single physics on several subdomains to the coupling of different sets of equations, like the case of fluid-structure interactions. The only requirements are that the different sets of equations are solved with the same iterative solver for 2 the same variable (e.g. velocity), as the coupling is carried out at the SpMV level. Figure 1: Example of applications of the proposed methodology. The parallel version of a SpMV is constructed in such a way to give the same result as its sequential counterpart. The way parallel SpMV is carried out depends on whether full-row or partial-row matrices are built on each MPI partition [6]. On the one hand, if full row matrices are considered, the nodes and associated matrix rows are assigned exclusively to one MPI partition. To construct such matrices, halo elements or extra-communications are thus required. In the finite element context, partial-row matrices are quite common as local matrices are assembled from element matrices coming from a partition into disjoint sets of elements. The rows of the duplicated interface nodes are thus not fully assembled on the interfaces. In this case, parallelization of the SpMV consists in exchanging the local results of the SpMV to assemble the interface contributions coming from the neighbors. This work extends this technique to cases where the nodes on the interfaces do not coincide. The proposed methodology thus generalizes the concept of coupling, for matching and non-matching meshes. 3 To introduce the method, we will start by considering the particular case of two subdomains. In Section 2, we will present the classical parallelization strategy used when dealing with partial-row local matrices and matching meshes. This strategy will be reinterpreted in terms of a Domain Decomposition method (DD), based on a Dirichlet/Neumann coupling to couple the different local meshes. In Section 3, this DD framework will enable us to devise the same strategy for non-matching meshes, by introducing transmission matrices to transfer Dirichlet and Neumann data from one mesh to the other. The resulting methodology can thus be used for both couplings between matching meshes, where each mesh represent a partition, and between non-matching meshes. In this context, the SpMV for matching meshes is just a particular case of this extended SpMV. Then, we will give in Section 5 some hints on how other operations of iterative solvers can integrate the nonmatching context as well. We will finally generalize the proposed strategy to an arbitrary number of subdomains in Section 4, by integrating in the same framework matching subdomains (as in the case of classical parallelization techniques) and non-matching subdomains as proposed here. The method will be illustrated through the analysis of the trace of a typical SpMV. The algorithms presented in this work will be given in such a way that operations can be carried out locally on each subdomain and when needed, communications will be indicated. The proposed framework should thus enable one to introduce non-matching meshes coupling in an already existing parallel code, where point-to-point communications are the main glue to obtain a global solution. 4 2. Revisiting the parallel SpMV In this section we will reinterpret the parallel SpMV in the context of domain decomposition methods. This framework will help us to devise the version for non-matching degrees of freedom coming from the gluing of nonconforming meshes. This will be done in next section. First of all, we will describe the classical strategy to parallelize the matrix-vector product. For this, we consider the following algebraic system coming from the discretization of a Partial Differential Equation (PDE) in a domain Ω: Au =b.(1) 2.1. Classical SpMV For the sake of simplification, we consider a partitioning of Ω into two subdomains Ω1and Ω2.Wedenotebyu1and u2the vectors of interior unknowns of Ω1and Ω2respectively, excluding the interface vector of unknowns denoted by uΓ. By performing a simple node reordering, System (1) can be written as:      A11 0A 1Γ 0A 22 A2Γ AΓ1 AΓ2 AΓΓ           u1 u2 uΓ      =     b1 b2 bΓ      .(2) In a parallel context, by considering a partition into disjoint subsets of elements, interface nodes are duplicated and the latter system is never fully assembled. See Figure 2 (Left). Instead, we obtain two independent systems 5 Figure 2: Two-subdomain decomposition with matching and non-matching meshes on the interface. for each subdomains:   A11 A1Γ AΓ1 A(1) ΓΓ    u1 u(1) Γ =  b1 b(1) Γ ,   A22 A2Γ AΓ2 A(2) ΓΓ    u2 u(2) Γ =  b1 b(2) Γ , where we have AΓΓ =A(1) ΓΓ +A(2) ΓΓ and bΓ=b(1) Γ+b(2) Γ.Fromnowon,we will append a superscript (i) to indicate a partial matrix or vector obtained locally in subdomain i. Let us consider the matrix-vector product y=Ax in this two-domain context. The parallelization of this product is based on the distributivity property of the multiplication. On the one hand, the product for interior 6 nodes is straightforward as each subdomain is able to perform it independently, that is yi=Aiixi+AiΓx(i) Γfor i=1,2, provided x(1) Γ=x(2) Γ=xΓ. For the interface we have: yΓ=AΓ1x1+AΓ2x2+AΓΓxΓ, =(AΓ1x1+A(2) ΓΓxΓ)+(AΓ2x2+A(2) ΓΓxΓ), =y(1) Γ+y(2) Γ,(3) where y(i) Γcan be calculated independently for i=1,2. In the MPI context, the assembled result on the interface yΓis obtained by exchanging the local values y(1) Γbetween neighbors, through the MPI non-blocking functions MPI Isend and MPI Irecv. A classical and non-optimized implementation is shown in Algorithm 1. In practice, communication of the interface values Algorithm 1 SpMV for matching meshes. 1: Compute local results for i=1,2: yi=Aiixi+AiΓxΓ, y(i) Γ=AΓixi+A(i) ΓΓxΓ. 2: Exchange results y(1) Γand y(2) Γbetween subdomains. 3: Assemble results in subdomains 1 and 2: yΓ=y(1) Γ+y(2) Γ. y(i) Γcan be overlapped with the computation of the interior result yi[7]. 7 2.2. SpMV as a domain (de)composition method The parallelization technique used previously to obtain the right answer on the interface can be reinterpreted in terms of a domain decomposition method. Let us duplicate formally the interface Γ into Γ1and Γ2and its associated unknowns, uΓ1and uΓ2, respectively. We then set the following Dirichlet/Neumann problem at the algebraic level [5]:   A11 A1Γ1 0I    u1 uΓ1 =  b1 uΓ2    A22 A2Γ2 AΓ22AΓ2Γ2   u2 uΓ2 =  b2 bΓ2 +  0 rΓ1  (4) where the residual is rΓ1=bΓ1−AΓ11u1−AΓ1Γ1uΓ1. The Dirichlet condition is the second equation that states that uΓ1=uΓ2 and the Neumann condition consists in assembling Neumann data rΓ1of subdomain 1 in subdomain 2 interface equation. We can easily check that this system is equivalent to system (2). Last system can be rewritten as:         A11 A1Γ100 0I0−I 00A 22 A2Γ2 AΓ11AΓ1Γ1AΓ22AΓ2Γ2                 u1 uΓ1 u2 uΓ2         =         b1 0 b2 bΓ2+bΓ1         A matrix vector y=Ax product using this matrix gives on the interface: yΓ2=(AΓ11x1+AΓ1Γ1xΓ1)+(AΓ22x2+AΓ2Γ2xΓ2), yΓ1=yΓ2. 8 In practice, this SpMV can be carried out exactly as in Equation 3, by doing local matrix-vector products and then exchange the solution. We have briefly shown how the substructuring method classically used to implement the parallel SpMV can be reinterpreted in terms of a Dirichlet/Neumann method at the algebraic level. Based on this domain decomposition framework, we devise in next section an equivalent method for nonmatching meshes, by introducing transmission matrices. Finally, it should be stressed that despite the fact that we use the formalism of the Dirichlet/Neumann method, the proposed implementation has nothingtodowiththeclassicalimplementations presented in the litterature. In general, the method is implemented in a staggered way (alaJacobi o a la Gauss-Seidel), sometimes as a preconditioner, sometimes as solver [8]. In our case the method is implicit and the associated matrix is exactly the same as the monolithic one, upon elimination of uΓ1. Thus the convergence of any iterative solver will be the same as the monolithic case, contrary to the case of staggered methods where the coupling does not even necessarily converge. 3. Extended SpMV to non-matching meshes We now consider non-matching meshes on the interface, as depicted in the right part of Figure 2. For the sake of calrity, we still consider only two subdomains. 3.1. Domain decomposition framework As degrees of freedom do not coincide, we need to introduce some transmission matrices to express the couplings between the unknowns on the subdomain interfaces. Let us introduce transmission matrices for the Dirichlet 9 In general, Krylov solvers [9] are based on simple solution updates like: uk+1 =uk+αpk. This equation corresponds to the solution update of the classical Conjugate Gradient method. Let us consider the first iteration with k=1: p0=b−Au0, u1=u0+αp0, and compute the updates on the interfaces. We need two initial conditions which satisfy the following relations to make it work: 1. RHS satisfies the Dirichlet condition: bΓ1 =TDbΓ2, 2. Initial solution satisfies the Dirichlet condition: u0 Γ1 =TDu0 Γ2.(10) We have p0 Γ1 =bΓ1−Au0|Γ1, =TDbΓ2−TD(Au0|Γ2), =TD(bΓ2−Au0|Γ2), =TDp0 Γ2,(11) where we have used Equation (10) and step 5 of Algorithm 2, or, equivalently, step 3 of Algorithm 3. We observe that the Dirichlet condition is directly transferred to the conjugate direction p. Then the update of the solution u 16 on Γ1yields: u1 Γ1 =u0 Γ1 +αp0 Γ1, =TDu0 Γ2 +αTDp0 Γ2, =TD(u0 Γ2 +αp0 Γ2), =TDu1 Γ2, where we have used Equations (10) and (11). Thus, the first update satisfies the Dirichlet condition as well. We could show that, recursively, the new updates will always satisfy the Dirichlet condition. Key message. If meshes are coinciding on the interface, then TD=TN= I, and we recover the classical parallel implementation of the SpMV given by Equation (3) with yΓ1=yΓ2=yΓ. The proposed technique is thus a generalization of the parallel SpMV for non-matching meshes. 4. Extension to an arbitrary number of subdomains Let us consider a parallel context with nMPI MPI partitions, referred to as parallelization subdomains, and which consist of disjoint sets of elements. Let us consider also ncou coupling subdomains, possibly non-matching, and coupled through Dirichlet and Neumann conditions. The intersections between the parallelization and coupling subdomains define a more general disjoint set of elements, simply referred to as subdomain. Figure 5 shows an example of four parallelization subdomains involving a non-matching coupling between two coupling subdomains. In practice, we start with an MPI partition into nMPI parallelization subdomains. Then, each MPI partitions is divided into 17 Figure 5: Example of four parallelization subdomain and two coupling subdomains. subdomains if it involves elements belonging to several coupling subdomains. In the example of the figure, we end up with six subdomains. To generalize the coupling algorithm for parallelization and coupling subdomains, we have to discriminate between Dirichlet and Neumann nodes. In this general context, the relation between Neumann and Dirichlet nodes should be understood as a master-worker relation, the master being the Neumann node, as shown in Algorithm 2: Neumann nodes accumulate residuals, and Dirichlet nodes transmit these residuals. Figure 6 shows the process for selecting Dirichlet and Neumann nodes when parallelization and coupling subdomains coexist. 18 Figure 6: Selection of Dirichlet and Neumann nodes in two steps. Step 1: selection of parallelization Dirichet and Neumann nodes. Let us start with the couplings between the parallelization subdomains, used in a classical parallelization strategy. The Dirichlet and Neumann nodes discrimination consists in selecting one of the duplicated nodes on the interfaces, mark it as Neumann, and then mark the duplicates as Dirichlet nodes. This is an explicit or implicit common practice in parallel codes, as illustrated on the left part of the figure. In fact, for the SpMV, the relation between duplicated interface nodes does not need to be explicit, as transmission matrices are equal to identity (see Algorithm 1). The distinction between these nodes is only necessary in the scalar product, where only the contributions of the interior and Neumann nodes are required to avoid duplication of the scalar product on the interfaces, as explained in Section 5.1. 19 Step 2: selection of coupling Dirichet and Neumann nodes. To fully establish a coupling for non-matching meshes, we have to select the Dirichlet and Neumann sides as well as the algorithms to compute the transmission matrices (interpolation, projection, transpose). The second step of the Dirichlet and Neumann node selection consists in marking as Dirichlet nodes the nodes located on the Dirichlet side of the coupling interface. In particular, Neumann nodes selected in the first step and located on the Dirichlet side must be converted into Dirichlet nodes (see the bottom marked zone on the right part of Figure 6). On the Neumann side, we then mark all the non-Dirichlet nodes as Neumann nodes (see the top marked zone on the right part of Figure 6). By performing steps 1 and 2, we thus ensure that Neumann nodes are uniquely defined on parallelization and coupling subdomains interfaces. The complete algorithm to couple the subdomains is very simple to describe, and a bit cumbersome to write formally. In brief, the Neumann nodes accumulate the residuals coming from its neighbors, through transmission matrices. The Dirichlet nodes are related to their Neumann counterparts through matrix TDand thus depend on the same neighbors as its master (Neumann) node. Nos let us go to the formal algorithm. Let nibe the number of neighbors of subdomain i, including itself, and nΓibe the size (number of nodes) of the interface of i. For the sake of simplicity, the transmission matrix TN ij which transmits data from subdomain jto subdomain iis dimensioned as nΓi×nΓj. In fact, only a sub-block of TN ij would be necessary as only subsets of Γiand Γjgive non-zero coefficients on 20 the matrix. The rows corresponding to Dirichlet nodes and the rows which do not involve neighbors jare set to zero. To account for its local residual, we also set TN ii =I. Using this formalism, the Neumann condition reads: yN Γi= ni  j=1 TN ij y(j) Γj.(12) As far as the Dirichlet transmission matrix TD ij is concerned, the rows corresponding to the Neumann nodes are set to zero, and we set TD ii =0. The Dirichlet condition on Γicoming from subdomain jis thus yD Γij =TD ij yN Γj=TD ij nj  k=1 TN jky(k) Γk, where we have used Equation (12). Now we need to sum over all the neighbors of iin order to take into account all the Dirichlet conditions coming from all the Dirichlet neighbors so that yD Γi= ni  j=1 yD Γij = ni  j=1 TD ij nj  k=1 TN jky(k) Γk. Therefore, as TN ij has no null rows on Neumann nodes and TD ij has no null rows on Dirichlet nodes, we can sum up yN Γiand yD Γito obtain yΓi= yN Γi+yD Γi. We finally end up with algorithm 4 which describes the extended SpMV for parallelization and coupling subdomains, and which consists of a generalization of Algorithm 3. This algorithm illustrates the extended SpMV, but this should not be implemented as is. In fact, when receiving the contribution of a neighbor j(step 3), not all the interface result of jis needed, but only that where columns of the transmission matrices are non-zero. In addition, just like in 21 Algorithm 4 Extended SpMV. 1: Compute local results on each MPI partition i: y(i) i=Aiixi+AiΓixΓi, y(i) Γi=AΓiixi+AΓiΓixΓi. 2: Send y(i) Γito all neighbors j. 3: Receive y(j) Γjfrom all neighbors j. 4: Assemble the result on the interface Γi: yΓi= ni  j=1 TN ij y(j) Γji + ni  j=1 TD ij nj  k=1 TN jky(k) Γk. classical parallelization techniques, interface communications represented in steps 2 and 3 can be overlapped with the SpMV for interior nodes using the non-blocking MPI communication subroutines. Let us analyze the performance of Algorithm 4. We consider the geometry illustrated in Figure 7 and compare the SpMV for three cases, all on 64 CPUs. Case 1 consists of a one domain simulation with only parallelization coupling. Case 2 involves two separate domains with only parallelization coupling. Finally, case 3 involves the same two separate domains, with parallelization coupling in each, and non-matching coupling between then. For the three cases, the total number of elements and nodes is around 4M, which results in an average of 62500 nodes per MPI partitions. By comparing case 3 with 2, we wish to measure the extra-communications due to the coupling. By comparing case 3 with case 1, we wish to point out the fact that the 22 Figure 7: Example to test the performance of the extended SpMV. (Top) Geometries and mesh sizes. (Mid.) Connectivities between the different MPI partitions due to parallelization and coupling. (Bot.) Partitioning. 23 partitioner (herein METIS[10]) sees two non-connected geometries and is thus unable to minimize the interface sizes on the non-matching meshes. See for example [11] for coupling aware partitioning techniques. On the middle part of the figure, we show the connectivities between the different subdomains. Each sphere represents a subdomain and its location corresponds to the subdomain center of gravity, while the lines symbolize connections between two subdomains. We can observe the extra connections between the MPI partitions in case 3 with respect to case 2. By taking a look at the MPI partititions, we also observe that the partitions of cases 2 and 3 are not aware of the coupling. In fact, in this work, the transmission matrices are computed in parallel after the partitioning. Figure 8 shows three different traces of one single SpMV, representing the different tasks (colors) carried out by the different preocesses (y-axis) along time (x-axis). The blue color represents the SpMV on the interfaces (Dirichlet and Neumann). The white color represents the non-blocking communications using the MPI functions MPI ISend and MPI IRecv. The red color, which dominates the computation, is the SpMV for the interior nodes. Finally, the purple represents the MPI Waitall which finalizes the SpMV. The top trace is that of the one-domain problem (case 1), the middle one the twodomain problem without coupling (case 2) and the bottom one the twodomain problem with coupling (case 3). We can observe the effects of the extra-communications involved in case 3. These communications do not exist in case 2, but they do in case 1. However, in this last case, METIS was able to take into account the minimization of the interfaces and thus to reduce 24 Figure 8: Traces of the SpMV on the geometries presented in Figure 7. The size of the window represents 2.7 ms. From Top to Bot.: Case 1: uncoupled problem represented by Algorithm 1; Case 2: coupled problem using a two-step communication scheme represented by Algorithm 2; Case 3: coupled problem using the one-step communication scheme represented by Algorithm 4. 25 286. doi:10.1007/978-94-011-1810-1 17. URL https://doi.org/10.1007/978-94-011-1810-1 17 [2] C. Bernardi, Y. Maday, F. Rapetti, Basics and some applications of the mortar element method, GAMM-Mitt 28 (2) (2005) 97–123. doi:10.1002/gamm.201490020. URL http://dx.doi.org/10.1002/gamm.201490020 [3] C. Farhat, F.-X. Roux, A method of finite element tearing and interconnecting and its parallel solution algorithm, Int. J. Num. Meth. Eng. 32 (6) (1991) 1205–1227. doi:10.1002/nme.1620320604. URL http://dx.doi.org/10.1002/nme.1620320604 [4] C. Lacour, Y. Maday, Two different approaches for matching nonconforming grids: the mortar element method and the feti method, BIT Numer. Math. 37 (3) (1997) 720–739. doi:10.1007/BF02510249. URL https://link.springer.com/article/10.1007/BF02510249 [5] G. Houzeaux, J. Cajas, M. Discacciati, B. Eguzkitza, A. GargalloPeir´o, M. Rivero, M. V´azquez, Domain decomposition methods for domain composition purpose: Chimera, overset, gluing and sliding mesh methods, Arch. Comp. Meth. Eng. 24 (4) (2017) 1033–1070. doi:10.1007/s11831-016-9198-8. URL https://link.springer.com/article/10.1007/s11831-016-9193-0 [6] G. Houzeaux, R. Borrell, Y. Fournier, M-Garcia-Gasulla, J. G¨obbert, E.Hachem,V.Mehta,Y.Mesri,H.Owen,M.V´azquez, Computational 32 Fluid Dynamics, Adela Ionescu Edition, Intech, 2017, Ch. High performance computing: dos and don’ts. [7] F. Magoul`es, F.-X. Roux, G. Houzeaux, Parallel Scientific Computing, Computer Engineering Series, Wiley-ISTE, 2015. URL http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1848215819.html [8] Domain decomposition methods for partial differential equations, Numerical mathematics and scientific computation, Oxford University Press, New York, 1999. [9] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, 2003. doi:10.1137/1.9780898718003. URL https://doi.org/10.1137/1.9780898718003 [10] K. Lab, Metis - serial graph partitioning and fill-reducing matrix ordering. URL http://glaros.dtc.umn.edu/gkhome/metis/metis/overview [11] M. Predari, Load Balancing for Parallel Coupled Simulations, Theses, Universit´e de Bordeaux, LaBRI ; Inria Bordeaux Sud-Ouest (Dec. 2016). URL https://hal.inria.fr/tel-01518956 33