# Multi-GPU Acceleration for Finite Element Analysis in Structural Mechanics
## 1 Introduction
We usually should avoid direct solvers due to the memory requirements for large systems of equations. In such cases, Krylov subspace iteration methods save memory consumption at the cost of increasing the computational burden, which can be alleviated using distributed computing [19; 20]. A key point to achieve high-performance GPU computing using Krylov methods [21; 22; 23] is the efficient implementation of the SpMV multiplication using sparse-matrix representation [24; 25; 26; 27]. We can obtain significant computing performance improvements by exploiting the fine-grained parallelism of the many-core architecture [28] and considering limited transfer rates over data channels, concurrency, and coalesced memory access. Another issue is the assembly process, which is highly parallelizable and can achieve significant speedups using GPU computing [29; 30], but it can require a high amount of memory. For this reason, assembly-free methods have been revisited recently for GPU computing [31] in finite element analysis [32; 33] and topology optimization [34] problems by exploiting the data locality properly. However, this strategy is inappropriate for finite element models using unstructured meshes, where we achieve poor or null speedup [35] using GPU computing. There is consensus that multigrid methods are the most efficient and popular techniques for solving large linear systems of equations. They exploit the smoothing property and the coarse grid principle. The former reduces the high-frequency error components, whereas the latter approximates the low-frequency error components on coarser grids, which we then prolongate to the finer grids. We usually classify these solving methods as geometric multigrid (GMG) and algebraic multigrid (AMG) approaches. The former uses the geometry and physics of the problem, whereas the latter only requires the coefficient matrix of the linear system of equations resulting from the discretization of Partial Differential Equations (PDEs). Regular grids facilitate the efficient calculation of the transfer operators between grids in GMG methods, which can be unaffordable in complex geometries. On the other hand, AMG methods are more flexible since they only use the information provided by the matrix of coefficients. For this reason, their use as a \"black-box\" function in finite element codes is prevalent. We also can classify AMG methods into two groups: classical AMG and aggregation AMG. The former obtains the hierarchical multigrid by partitioning the nodes into coarse and fine nodes and defines the interpolation operator (or transpose of restriction operator) by the weighted sum of such nodes [36]. In contrast, aggregation AMG methods [37] aggregate a few fine grid nodes to form a coarse grid node and define the interpolation as a piecewise constant operator, which results in an interpolation operator with only one zero per row, reducing the memory requirements meaningfully in comparison with classical AMG [38] and improving the efficiency of the interpolation operator. These features have a significant beneficial effect on the performance of GPU computing. Multigrid methods are commonly used as a preconditioner of Krylov subspace methods [39; 40] in the structural mechanics context. The underlying idea is that the interpolation operator of the multigrid method will hardly be optimal, which makes it less efficient for some specific error components. When this occurs, the convergence is slow despite almost all the error components reducing quickly. Krylov subspace methods usually eliminate these error components efficiently [41]. Commonly, this is a more efficient solution than improving the construction of the interpolation operator. In this context, GPU computing using geometric multigrid methods (GMGs) as a preconditioner has shown clear benefits, both in structural analysis [42] and topology optimization [43; 44; 45].