GPU acceleration of equations assembly in finite elements method-preliminary results

Abstract

The finite element method (FEM) is widely used for numerical solution of partial differential equations. Two computationally expensive tasks have to be performed in FEM-equations assembly and solution of the system of equations. We present mapping of the equations assembly problem for StVenant-Kirchhoff material to GPU computation model and show results of its early implementations outperforming our single core CPU implementation by factor of 15. Moreover, presented method used for equations assembly problem solution is more general and can be considered as a technique to manage complex composition of medium-grained operations with different requirements of GPU resources. I. EQUATIONS ASSEMBLY In this paper, CUDA-based acceleration of equations assembly in static deformation FE model is presented. The model coupling the applied loads f and resulting deformation u is based on the theory of elasticity. Beside the static equilibrium relating the internal stress tensor σ with the applied loads f , the model is given by the formulation of strain tensor. In our case, non-linear Green strain tensor γ is used: γ ij = 1 2 (∂ j u i + ∂ i u j + ∂ i u m ∂ j u m) (1) where ∂ i denotes derivative of the argument w. r. t. the i-th component. In our case, hyperelastic model is considered given by a stored energy function W such that σ ij = ∂W ∂γ ij. (2) For the StVenant-Kirchhoff material, the stored energy W is defined as W = λ 2 γ ii γ jj + µγ ij γ ji (3) where λ and µ are material coefficients known as Lamé constants. The resulting boundary value problem (BVP) is given by a governing equation which is built from the equations given above defined over a domain Ω and boundary conditions. The BVP cannot be solved analytically due to the non-linearity of the governing equation and irregularity of the domain. Therefore, finite element method is used to reformulate the problem using weak formulation over the domain which is discretized into simple-shaped elements such as cubes and tetrahedra. More precise formulation of FEM is beyond the scope of this paper (see e. g. [1]). The final outcome of the method is represented by a non-linear system K(u) = f of algebraic equations, which is usually solved by some iterative method such as Newton-Raphson, where in each iteration, a linearized problem A(u)∆u = f −K(u) is solved using some actual estimation of u. In this paper, the procedure assembling the matrices A(u) and K(u) is studied. It is supposed that the domain is discretized by linear tetrahedral elements having four vertices each. During the assembly procedure, small 12 × 12 matrices A e (u) and K e (u) are assembled for each element e. Denoting φ the linear shape functions provided by the finite element method and δ ij the Kronecker δ-symbol, the matrices A e (u) and K e (u) are defined by components using the tensor notation with summation convention as follows: A N x = Z Ω e ∂W ∂γij (δxi + ∂iux) ∂N φjdV − Z Ω e fxφN dV − Z ∂Ω e gxφN dS (4) K N M xy = Z Ω e  ∂ 2 W ∂γij∂γ kl [(δyi + ∂iuy) ∂jφN ] [(δ xk + ∂ k ux) ∂ l φM ] + δxy ∂W ∂γij ∂iφN ∂jφM ff dx. (5) where the indices M, N = 1. .. 4 run over the vertices of the element and x, y = 1. .. 3 are the components of the 3D problem. The integration over the volume Ω or its boundary ∂Ω is performed numerically using Gauss quadrature. Once the system of equations is assembled, it has to be solved. We focus on equations assembly problem only, as solving of both sparse [2] and dense [3] systems of equations has been successfully accelerated using GPUs. The additional advantage of performing equations assembly in GPU is that the assembled system is already in GPU's memory eliminating the PCI-E overhead. II. PROBLEM MAPPING TO GPU COMPUTATION MODEL In this chapter, the mapping of equations assembly problem to GPU computation model is discused. It is out of the scope of this paper to describe CUDA-enabled GPU programming and memory model, therefore we refer to NVIDIA CUDA Programming Guide [4]. The system of the equations solved in FEM is composed from results of equations assembly for all elements of the mesh. The amount of computations needed for the equations assembly for one element is insufficient to utilize GPU well,