Numerical integration on GPUs
for higher order finite elements

We consider the problem of numerical integration for higher order 3D finite element approximations. The main feature of the problem that we address is the large variation of required resources, processing time and memory, for different orders of approximation, i.e. different degrees of polynomials employed as finite element shape functions.

We investigate the characteristics of the problem and try to design a portable kernel for different problems and different processor architectures. Hence, we restrict ourselves to OpenCL programming model. It allows us to construct kernels for graphics processors supporting the model (currently NVIDIA, AMD and Intel hardware), as well as for CPUs. Our previous investigations [34] as well as results reported in other papers (e.g. [27]) show that OpenCL implementations achieve performance only several percent lower than native CUDA model for NVIDIA GPUs.

We situate our research in the context of creating a general purpose finite element framework ModFEM [35], for continuous and discontinuous Galerkin approximations. The code has modular structure in which different software components are responsible for different phases of finite element calculations [36]. Each component can exist in several variants and final applications for different problems and solution strategies are created using suitable compilation options.

The code has also a layered structure with separate modules responsible for shared memory computations, that form a core layer and that are further extended to distributed memory environments by suitable overlays [37, 38], belonging to a different layer. The code operates on unstructured adaptive meshes with proper load balancing for parallel execution [39]. The code has abilities to run on clusters with message passing and is currently ported to different new computer architectures, such as GPUs and heterogeneous processors, like CBE [40].

In each particular code created by ModFEM framework there is a separate approximation module that is responsible for the creation of element stiffness matrices. The matrices are further transferred to linear equations solvers using a generic interface [41]. Thanks to the layered structure of the code there is only one numerical integration routine, in the core layer, for both, shared and distributed memory environments. The interface between numerical integration procedure and the solver of linear equations allows for employing different, direct and iterative, solvers of linear equations, possibly using thin adapters [42].

The same integration routine is also used for different types of meshes and problems solved. This approach proved successful in many application domains [43, 44, 45]. As a main problem of the current investigations we pose the problem whether it possible to maintain this portability and genericity by creating a single kernel executed on GPUs for different problems, meshes and computer architectures.

We consider numerical integration only, leaving the problem of assembly (and also application of boundary conditions) to different software and possibly hardware components. In our view this smaller grain of considered GPU code gives more flexibility in designing solution strategies, where CPU and GPU cores can work together, each group performing the tasks it is most suitable for [17]. In order to maximize this flexibility of the approach we consider different options for passing the results of calculations. The algorithms presented in the paper do not prevent the use of different assembly techniques, such as e.g. presented in the discussed above studies devoted to the subject.

We attempt to investigate the possibility of creating a general purpose numerical integration procedure executed on graphics processors that can either replace or complement the CPU integration routine. We do not attempt to create GPU-only code. Considering relatively small sizes of memory for current accelerators with GPUs, in the order of several GBytes, this would limit the sizes of problems that we solve. Since we aim at large scale problems we want to use large memories, in the range of hundreds of GBytes, for each shared memory node.

The use of GPUs as accelerators for CPUs, corresponds also to our distributed memory strategy where we assume two separate levels of decomposition. GPU calculations fit into shared memory level calculations controlled by CPU code and due to the design philosophy can be directly used in distributed memory calculations, based on CPU code and its message passing overlays.

Following the arguments presented above, we assume that our general purpose finite element numerical integration procedure does not require double precision calculations and we restrict our interest to single precision implementations. We plan to consider extensions to double precision for particular application domains that would pose such requirements.

The investigations presented in the paper form the continuation of an earlier work. In [52] we presented numerical integration on GPUs for 2D linear elasticity, CUDA environment and older generations of GPUs. This work was extended to OpenCL programming model in [34] and to 3D Laplace operators in [53]. We investigate separately numerical integration procedures for CBE processors in [52] and [54].

The paper is organised in the following way:

• 1.

  Section 2 defines the problem of numerical integration for finite element approximations

• 2.

  Section 3 describes the programming model that we use and its mapping to selected GPU architectures

• 3.

  Section 4 contains details of several versions of numerical integration procedures for GPUs that we developed

• 4.

  Section 5 presents results of tests performed to assess the performance of the code

• 5.

  Section 6 contains final remarks.

CUDA [48] and OpenCL [49, 50] are the two most popular programming models for GPUs. CUDA was the first platform for programming GPUs that gained widespread acceptance and popularity. Still it is the most advanced programming environment, aimed at getting the most of graphics hardware produced by NVIDIA.

The OpenCL programming model and environment was created later, however, from the beginning it targeted a broader spectrum of hardware, including, apart from GPUs, standard multi-core CPUs and hybrid architectures, like Cell Broadband Engine (CBE). In the context of creating a general purpose, portable code for finite element simulations we concentrate on the OpenCL programming model.

Despite different targeted architectures, CUDA and OpenCL are very similar. Most of CUDA codes, that do not use advanced capabilities of the computing environment, can be easily transferred to OpenCL. Most of notions used in one programming model have direct analogues in the other. Software design presented in a sufficiently generic manner for one of the models can be considered as a design for both. We believe that the procedures that we design for OpenCL can be directly ported to CUDA, with no performance loss. They can be further optimized for NVIDIA hardware using techniques specific to CUDA.

In order to make our presentation easily accessible to persons more accustomed to CUDA programming, when introducing OpenCL notions, we recall in parantheses respective CUDA names if they differ from the ones used in OpenCL. In the description of GPU algorithms we use a mixture of the two styles, that we consider intuitively clear.

OpenCL code is compiled and run on a given platform, that represents the environment for code execution. Each platform provides a set of devices (CPUs, GPUs, accelerators). For each platform there exists a host system that runs standard code and manages the execution of a set of kernels on devices. Kernels are functions written in the OpenCL language, a slightly modified version of the C language. Since the specification does not provide separate notions, by kernels we will understand both, functions in the form of source code, as well as in binary form that is send for actual execution.

Devices are composed of compute units, which are further divided into processing elements. Individual threads of execution (called work-items in OpenCL specification – the notion that we will not use in the paper) are running on processing elements, that are in principle capable of executing scalar and vector instructions.

For running on compute units, threads are grouped into work-groups (thread-blocks in CUDA nomenclature), with one work-group assigned to a single compute unite (each compute unit can, however, have many work-groups assigned to it). Threads within a single work-group execute concurrently, can share some data in fast memory and can be synchronised using fast system calls. Different work-groups are scheduled independently.

For running on a device, thread work-groups are further combined into the whole set of threads executing a kernel. Threads within a work-group have local identifiers, that may have form of 1, 2 or 3 dimensional arrays. Global identifiers for threads within the whole executing set can have similar form. We use simple one dimensional local and global numbering of threads. Work-groups have also their indices, and, in the same way as in the case of threads, we use one dimensional numbering of work-groups.

OpenCL execution model specifies the events that has to occur in order for device code to be run. The first phase includes initialization of the specific OpenCL data structures and checking the available devices. Then the code to be executed on the devices has to be read from source (a distinguishing feature of OpenCL) or binary file, possibly compiled and prepared in the form of a kernel. During compilation several options can be passed to the OpenCL compiler.

Kernels can accept arguments, hence, before executing a kernel, the space for arguments has to be allocated on device and arguments send to the device memory. The host code can allocate space for variables and arrays in different explicitly available to programmers types of memory.

The memory transactions are performed using a typical for OpenCL strategy, by sending a request to the OpenCL management layer. The requests are then realized asynchronously to the host code. In the same way the requests for kernel execution and transferring back data from device memory to the host memory are realized.

OpenCL memory model defines several types of memory regions explicitly available to programmers. Different memory objects can be created for OpenCL kernels with different mappings to hardware resources. We briefly describe simple classification of variables that we use in our calculations.

Individual variables defined in kernels belong to private (local in the CUDA nomenclature) memory. Each thread has its own copy of each variable, preferably stored in registers.

Assignment to other memory regions is achieved through specific qualifiers. We use the arrays stored in the following types of memory:

• 1.

  global – such arrays are visible to all threads executing the kernel

• 2.

  constant – being a part of global memory with read only access for threads

• 3.

  local – arrays in this fast memory are shared by threads in a single work-group

Since the notion of local memory has different meanings in OpenCL and CUDA terminology, we do not use it in the rest of the paper. Instead, for OpenCL local memory we use the CUDA term ”shared memory” that in our view better reflects its character.

Since OpenCL is aimed at portability of created codes, it contains procedures that allows for adapting to different platforms and devices. The code can query the environment to get information on many available resources. We employ several of such functions, the ones that we use will be described in the section devoted to parallel GPU implementation of numerical integration.

In this section we briefly describe the hardware resources available to CUDA and OpenCL programmers. Since the GPU hardware changes rapidly, we do not give technical details, that are different for different architectures and different generations of GPUs. We try to provide a generic description that should serve as a justification for several of design choices that we make. We restrict the description to NVIDIA and AMD GPUs as the two most popular families.

Explicitly defined OpenCL and CUDA abstractions are mapped to different hardware resources for particular GPU architectures. Compute units are called streaming multiprocessors for NVIDIA hardware, processing elements are often called CUDA cores for NVIDIA and stream cores for AMD. The details of GPU organization can be complex and the mapping of compute units and processing elements to hardware units may differ for different architectures [48, 51].

The memory hierarchy is implemented in AMD and NVIDIA architectures in a similar way. There are registers, off-chip DRAM for global and constant memory and on-chip hardware for shared and cache memory (constant memory is cached, the global memory is cached in recent generations of GPUs). The access times and bandwidths for different memory regions vary for different architectures. The variations include also relative speeds of different memory types. Specifically, shared memory may be almost as fast as registers for certain architectures (e.g. NVIDIA Fermi) or several times slower (e.g. AMD Evergreen).

One important thing to mention, is handling of automatic variables, private to each thread. They are preferably stored in registers, but if too many registers are required by a kernel, either the kernel is not executed or register spilling occurs. For some architectures, spilled variables are stored in global memory, that significantly slows down the execution of a kernel. Recent GPUs utilize cache memory for variables spilled from registers, that makes the problem of minimizing the number of registers in a kernel less crucial.

We use the term GPU memory for the global, device memory of GPUs. For each architecture that we consider, GPU memory is separate from the host memory and connected to it through PCIe bus. Slow PCIe connection should be taken into account when designing GPU kernels, as well as the fact that memory initialization forms another source of overhead for GPU execution.

The most important fact related to mapping software to hardware and not explicitly visible to programmers is the way threads are scheduled for execution. Both NVIDIA and AMD use SIMD style, where threads are grouped into sets that are called warps for NVIDIA and wavefronts for AMD. All threads in such a warp/wavefront execute one instruction at a time.

If the instruction executed by a warp/wavefront is a global memory access and the locations accessed by different threads do not form one of specific patterns, one memory access instruction is executed using many memory transactions. However the proper arrangement of accesses allows hardware to, as it is called in CUDA nomenclature, coalesce accesses to different locations into a single memory transaction.

Still such a transaction takes many clock cycles. In order to prevent pipeline stalls several warps/wavefronts have to be executed concurrently on a compute unit. We try to ensure in our kernels that many warps/wavefronts are scheduled concurrently, as a key ingredient for achieving high performance on GPUs.

The problem of proper arrangement of memory accesses concerns not only the global GPU memory, but also the shared memory of compute units, that is implemented in hardware using many memory banks that have to be accessed in a way that eliminates bank conflicts. One simple rule, that should ensure the proper memory usage, is to make threads access subsequent elements of arrays storing single precision floating point numbers.

We adopt this rule in the design of our kernels. Hence, whenever we talk about accessing vectors and matrices in memory, either global or shared, we mean accessing 1D arrays with subsequent entries read or written by subsequent threads.

It can be seen from the description above that the size of warps/wavefronts is an important characteristic of hardware. The number of threads in a work-group should be a multiple of warp/wavefront size. However the details of execution for different GPUs make the picture less clear [48, 51]. We adopt, as a rule of thumb, based on the descriptions of utilized architectures, that the size of work-groups in our kernels forms a multiple of 64 (the size of two warps for NVIDIA and one wavefront for recent AMD GPUs).

Apart from standard optimizations that we perform, such as removing dependencies, reducing the strength of expressions, etc. and the guidelines presented in the previous section, we utilize several further rules related to specific character of GPUs.

One rule, is to maximize the number of threads executing a kernel. This is done not only due to the fact that GPUs have many processing elements, but also the fact that in order to hide memory access latency the number of threads ready to be executed should be a multiple of the number of processing elements. By which factor however, it depends on particular GPU architecture. We assume that at least several threads per processing element should be designated.

This can be achieved in two ways. One is to increase the size of a single work-group. By doing this the number of threads for a single compute unit is also increased. The second way is to increase the number of work-groups. Different work-groups can be scheduled to the same or to different compute units, hence the total number of threads for a GPU grows.

Which strategy to choose depends on several factors. One of such factors is the amount of memory resources required by threads. Increasing the number of threads, increase the required number of registers. Increasing the number of work-groups increase the required amount of shared memory. This may suggest as a good compromise to increase the number of threads within one work-group, up to the optimal maximum for a given compute unit, and then increase the number of work-groups to spread as many threads as possible across many compute units. The strategy has to be accompanied by a proper kernel design with minimal possible register and shared memory usage.

As we mentioned already, in designing our kernels we aim primarily at portability of the code across different problems considered and different processor architectures. Hence, we do not describe here and do not use specific optimizations targeted at particular architectures [59, 50]. This leaves the space for further improvements, that can yield better execution performance in particular circumstances.

The design of our parallel GPU code for higher order numerical integration is based on several assumptions:

• 1.

  single kernel invocation is performed for a large set of finite elements – this is necessary since kernel invocation incur substantial overhead that has to be amortised over as many kernel operations as possible

• 2.

  one element is assigned to one work-group of threads, however one work-group operates on a set of elements, forming a subset of all elements associated with kernel execution

• 3.

  for a single work-group and a single element parallelization is performed in the data-parallel manner, with each thread performing calculations for a set of entries in the element stiffness matrix

Several further design decisions have to be made then. How many threads should be in a work-group? How many entries of $\mbox{\boldmath$A$}^{e}$ should a single thread process? How arrays of shape functions should be accessed by threads? All these questions should be answered based on characteristics of the performed calculations in the light of great variations of the sizes of data structures associated with different orders of approximation.

One of such decisions concerns the use of values of shape functions. For sequential codes and low order approximations it may happen that computing the values of shape functions ”in flight” may present no overhead as compared to reading them from memory. However for GPU execution, the facts that computing shape functions and their derivatives has small degree of concurrency (some parts of calculations are inherently sequential) and that these values are the same for all elements processed by the kernel, make precomputing the values and storing them in arrays in memory the only solution that can lead to high performance. Since the sizes of such arrays, presented in Table 2, exclude the use of some faster memory, they are stored in global GPU memory.

The size of the element stiffness matrix, at least for the higher considered orders of approximation, is also too big to fit into fast GPU memory and $\mbox{\boldmath$A$}^{e}$ has to be kept in the slower global memory. It would be inefficient to perform summation of contributions from different integration points in global memory, hence it is performed using temporary variables in fast memory. We consider two choices: the values are stored either in registers (individual for each thread) or in shared memory (common to a work-group).

If threads are using registers for storing computed values of $\mbox{\boldmath$A$}^{e}$ than for a single loop over integration points as many blocks of $\mbox{\boldmath$A$}^{e}$ as there are threads in a work-group are computed. After summation, the values are written to global memory and in the next step the threads has to consider new entries to be summed up. However this step requires repeating calculations in lines from 6 to 9 in Algorithm 2, including shape function reads and calculations related to Jacobian matrices, many of which are redundant. The whole procedure reduces to introducing another loop in Algorithm 2. Before the loop over integration points, there must appear a loop over parts of stiffness matrix processed at once by all threads in a work-group.

The number of parts into which the whole element stiffness matrix $\mbox{\boldmath$A$}^{e}$ is divided can be simply computed based on the stiffness matrix size and the number of threads in a work-group. The number of parts is equivalent to the number of redundant repetitions of shape function reads and calculations related to Jacobian matrices. This redundant calculations form an overhead as compared to sequential code and can lead to lower parallel speed-up. We investigate two strategies to make this overhead as small as possible.

The first strategy is simple and consists in making the calculations in lines 6-9 in Algorithm 2 as fast as possible. This includes the option for reading, instead of computing, the necessary terms related to Jacobian matrix and its inverse (we will call them later Jacobian terms). Since the calculations are inherently sequential this single thread section significantly slows down kernel calculations. However, the option of sending data implies that the values of Jacobian terms for all elements processed by the kernel and for all integration points are computed by the host code and sent to global GPU memory. These additional calculations and memory transfers form another overhead, the influence of which we will assess in computational experiments.

We will test two versions of the numerical integration algorithm related to storing blocks of $\mbox{\boldmath$A$}^{e}$ in registers, the first, denoted by REG_JAC with assumption that Jacobian matrices are computed on GPU, the second denoted by REG_NOJAC assuming that the values of Jacobian matrices are sent to GPU and further read by threads. Both versions are depicted as Algorithm 3. Whenever the operation of reading is specified without further details, it means a coalasced read from global to shared memory.

The second strategy for reducing single thread region overhead, mentioned already as an option for parallelization, is to allow threads to compute more blocks for a single quadrature point, by storing them not in registers but in shared memory. This has the advantage of allowing for more explicit management of resources, since in OpenCL the size of shared memory can be queried by the host code.

In the first strategy, the shared memory may be not fully utilized (it contains problem dependent coefficients, element geometry data and shape functions with derivatives for a single integration point). In the second approach the size of the part of shared memory devoted to store the part of stiffness matrix processed by a single work-group of threads can be explicitly computed, assuming the full utilization of shared memory and taking into account the intended number of active work-groups for a single compute unit. As the drawback of this approach the use of shared memory, that may be slower than registers (even assuming properly arranged accesses by all threads) has to be considered.

The two versions of the algorithm that use shared memory for storing computed parts of element stiffness matrices, SHM_JAC with Jacobian terms computed on device and SHM_NOJAC with Jacobian matrices sent to GPU global memory are presented as Algorithm 4. It should be noted, that because of the large size of the element stiffness matrix and usually the small size of the compute unit shared memory, even here, the loop over parts of the stiffness matrix has to be performed.

In Algorithms 3 and 4 not all threads are engaged in all kernel operations. In order not to waist the global memory, that may be used by other kernels operating on computed stiffness matrices, there is no padding in global arrays used by kernels. Since the global data is always accessed in the kernels in a coalesced way, the number of threads performing global read/write operations is specified exactly in the kernels, with some threads remaining idle.

Figure 1 presents the kernel workflow with indicated synchronization barriers (using horizontal bars) and the division of work among threads. Different numbers of threads perform different access operations, hence the term ”subset of threads” denotes different numbers of threads for each read/write operation.

Running the kernels described above requires the existence of a CPU code that manages, on the host side, the execution of GPU functions. We present the host code in the context of large scale calculations, where CPU code may represent one of processes in distributed calculations and the number of elements for single CPU manager thread is large (for higher orders of approximation large may mean tens of thousands, for low order approximations the number may be in the range of millions).

In order to prepare for numerical integration, the host code classifies all assigned to it elements from the mesh or submesh (in distributed calculations) into sets of elements of the same order of approximation. When assembly procedure requires this, each set can be further divided into subsets of elements that do not contribute to the same entries in the created global matrices and vectors (e.g. by some colouring algorithm). Further procedures that we describe concern one of such subsets.

The host code queries the platform on which calculations are performed about the properties of the GPU device. These include:

• 1.

  the size of GPU global memory (and the maximal size for a single allocated array)

• 2.

  the size of shared (local in OpenCL nomenclature) compute unit memory

• 3.

  the size of constant memory

• 4.

  the maximal number of threads in a work-group i.e. the maximal work-group size

• 5.

  the maximal number of computing units

• 6.

  the maximal number of threads for kernel

After that the CPU code calculates the parameters send as arguments for the GPU kernel. Based on the number of blocks in the element stiffness matrix and the maximal work-group size, the number of threads in a single work-group is computed (as a multiple of 64 as we mentioned in Section 3.3). This number, the size of work-group, is maximized, without exceeding other GPU resources and going beyond the number of blocks in $\mbox{\boldmath$A$}^{e}$.

Then, the host code computes the number of parts into which $\mbox{\boldmath$A$}^{e}$ is divided. For SHM algorithms it is calculated in such a way that the full available compute unit shared memory is used. This gives the size of the part of stiffness matrix and the number of blocks processed by a single thread in a sequence for the part (it is equal 1 in REG algorithms). Based on this number and the size of work-group the number of parts of $\mbox{\boldmath$A$}^{e}$ is computed.

In the next preparatory step, the number of elements processed by one kernel invocation is calculated. It is maximized based on the available size for an array storing element stiffness matrices in GPU memory. It is also adjusted to the number of elements in a subset of elements with the same degree of approximation and possibly the same colour.

As the final calculation, based on the number of elements processed by the kernel, the number of elements assigned to a single work-group (associated with the loop in line 2 in Algorithms 3 and 4) is computed. In order to maximize the number of active work-groups, the number of elements per work-group is minimized. However, still must be large enough so the number of threads does not exceed the limits specified for the GPU device.

After computing parameters for kernel invocation, data forming kernel arguments are created. This include:

• 1.

  the set of calculated or retrieved from finite element data structures parameters of execution – sent to GPU in an array

• 2.

  quadrature data – for JAC versions of algorithms

• 3.

  geometry related data (nodal coordinates) for elements – for JAC variants

• 4.

  Jacobian terms (determinant of the Jacobian matrix and 9 entries of the inverse Jacobian matrix) for all elements and all quadrature points – for NOJAC variants

• 5.

  array of problem dependent coefficients for elements (possibly combined together in one data structure with geometry data)

• 6.

  the values of all shape functions at all integration points for a reference element (or elements) corresponding to all elements processed by the kernel

The prepared data are sent to the GPU memory in several large transfers. The fact of associating as many finite elements as possible with a single kernel invocation should give performance gains, since the overhead associated with data allocation on device and memory transfers is amortized over many elements.

Finally some arrays are allocated in the GPU memory. For SHM variants the space for a part of stiffness matrix is allocated in the shared memory. For all variants the workspace for shape functions is allocated in the shared memory and the space for the calculated element stiffness matrices is allocated in the global memory.

The integration kernel is invoked in a usual OpenCL style with the number of threads in a single work-group and the total number of threads for kernel as arguments.

Computed stiffness matrices either may be copied back from GPU to host memory and transferred to the assembly routine or may be left in GPU memory as an input for some different kernel performing assembly. The first alternative may correspond to situations where there are more elements on the computational node than can fit into GPU global memory, and hence, the GPU performs numerical integration for another portion of elements. For the second alternative, not only assembly can be performed on GPU but the whole solution to the system of linear equations or an update of global vector of unknowns. In such cases computed stiffness matrices are not copied back to the host memory at all.

As a model problem for numerical examples we consider a 3D linear elastostatics. The differential formulation of linear elastostatics for the vector of displacements $u$ = ${u_{1},...,u_{d}}$, $d=3$ is:

where $c_{ijkl}$ are in this case elasticity coefficients and summation convention for repeated indices is used. Particular forms of matrices $c$ depend upon materials considered. The number of space dimensions $d$ represents, in our case, also the number of equations in the system and the number of components of the unknown function $u$.

After discretizing the computational domain $\Omega$ into finite elements and using a proper basis functions $\hat{\phi}^{r}$, $r=1,..,N$, (the same functions are used for all components of $u$) the approximated displacements $\mbox{\boldmath$u$}^{h}$ are expressed as linear combinations:

where $U_{k}^{r}$ are the discrete unknowns, the coefficients of linear combinations of basis functions.

The resulting final weak statement of the problem that forms a basis for computational procedures is:

The domain integral in equation (3) is computed as the sum of element integrals. This corresponds to the fact that the final system of linear equations is assembled from element stiffness matrices, the fact that we adopted as the starting point for the formulation of the numerical integration problem. The analytical expression for a single entry of an element stiffness matrix, in the case of our model problem of linear elasticity is the following:

where $\frac{\partial{\phi}^{r}}{\partial\mbox{\boldmath$\xi$}}$,$\frac{\partial{\phi}^{s}}{\partial\mbox{\boldmath$\xi$}}$ are vectors of derivatives of respective element shape functions with respect to local, reference element coordinates, $w^{Q}_{I}$ denotes the weights associated with particular integration points for a particular quadrature and $\mbox{\boldmath$J$}_{\mbox{\boldmath$T$}_{e}}$ denotes the Jacobian matrix of transformation ${\mbox{\boldmath$T$}_{e}}$ from a reference element to the real element $\Omega_{e}$ and $\frac{\partial\mbox{\boldmath$\xi$}}{\partial x_{i}}$ forms a column of the inverse Jacobian matrix $\mbox{\boldmath$J$}_{\mbox{\boldmath$T$}_{e}^{-1}}$ corresponding to the inverse transformation ${\mbox{\boldmath$T$}_{e}^{-1}}$). This expression directly corresponds to kernel Algorithms 3 and 4.

For our model problem of linear elastostatics we assume that PDE coefficients, being material data, can be different for each element, but are constant over a single element. In practice, we use only two parameters: Young modulus and Poisson ratio. The particular form of the array of coefficients implies the operations performed for a single $3\times 3$ block of an element stiffness matrix. There are 63 floating point operations necessary in the case of standard updates done by sequential codes.

The same number of operations for a single block is also performed by the GPU kernels, however, the number of blocks involved is larger, due to the fact that the number of threads in a work-group is a multiple of 64 and is maximized in our kernels. Hence, the number of operations actually performed by hardware is larger for GPUs than for standard processors.

For execution on GPUs apart from floating point operations, there are some operations performed by each thread in the innermost loop of Algorithms 3 and 4 due to the necessity of computing, for each block, its location in the element stiffness matrix, based on the thread ID and its work-group ID. These calculations involve integer division operations, expensive for GPUs.

The created kernels have been tested for the full range of considered orders of approximation and several GPUs. For comparison, the same integrations have been performed on a server with x86 processors.

We present the results for two GPUs on graphics cards: NVIDIA GTX580 and AMD HD5870. The first GPU represents the Fermi architecture [56], has 16 streaming multiprocessors and 512 CUDA cores. The GPU memory size is 1.5 GB, with the maximal bandwidth 192.4 GB/s. The number of single precision floating point registers available for threads on a single compute unit is 32786 and the single precision theoretical peak performance is 1581 GFlops. The parameters for the AMD GPU, belonging to the Evergreen family [51], are: 20 compute units, 1600 processing elements, 1 GB global memory, 153.6 GB/s memory bandwidth, 16384 registers and 2.72 TFlops single precision theoretical peak performance. For NVIDIA card we used CUDA 5.0 SDK that provided OpenCL library and for AMD we used APP (version 923.1) environment. Tests for both GPUs were done with Linux operating system and proper driver software provided by NVIDIA and AMD.

The x86 platform consisted of an eight-core Intel XEON E5-2670 processor (2.6 GHz, 332,8 double precision GFLops, 51.2 GB/s maximal memory bandwidth), Linux operating system and Intel C compiler. For comparison with GPUs we have used the standard version of the ModFEM code with double precision calculations and parallelization of numerical integration achieved by parallelizing the loop over elements using OpenMP.

It has to be noted that for comparison we have used an optimized version of CPU code (although for double precision calculations) and we employed a recently (2012) introduced, powerful x86 processor. The two GPU processors are relatively older, by two (NVIDIA) or three (AMD) years, which in the light of the Moore’s law and implied speed-ups, puts GPUs in less favourable position.

Kernel tests were performed for four variants of numerical integration algorithm. We have not introduced specific, hardware oriented optimizations. The same kernels were executed for each order of approximation and each GPU.

For each run the execution begins with retrieving data concerning the platform and the device utilized. These include the size of global memory and the number of compute units (for the values appearing in test runs, see hardware description section) as well as other characteristics, such as the maximal size for memory object allocation (384 MB for NVIDIA and 128 MB for AMD), the size of shared (local in OpenCL terminology) memory (48 kB for NVIDIA, 32 kB for AMD), constant memory size (64 kB for both devices) and the maximal number of threads (work-units) in a work-group (1024 for NVIDIA and 256 for AMD).

Based on the retrieved data, parameters of kernel execution are computed. These parameters are presented in Tables 3 and 4 for the considered orders of approximation $p$ from 2 to 7 and the two GPUs. The only important parameter not included in the table is the number of threads in each work-group. Since we aim at maximizing the number of threads running on a GPU, we use 256 threads per work-group for AMD and 512 threads per work-group for NVIDIA (except the case of $p=2$ where we use 192 threads because of small size of stiffness matrices). The fact that we do not use maximal number of threads per work-group for NVIDIA, results from the large register consumption of our kernels (in the range of 30-37 for AMD and 53-63 for NVIDIA), due to performing several nested loops and complex calculations for a single block of entries of $\mbox{\boldmath$A$}^{e}$.

We present performance results for our numerical integration kernels in terms of execution times and GFlops rates. The results in subsequent Tables 5, 6 and 7 show for subsequent processors, Xeon E5-2670, GeForce GTX580 and Radeon HD5870, the execution times for different phases of numerical integration procedure and the computed performance of calculations. Several remarks below explain the assumptions on which measurements and performance calculations are based.

• 1.

  Each presented table contains timing results for subsequent stages of the whole procedure of numerical integration. This reflects the possibility of arranging the stages in different ways. The initial stage is important when the first group of elements is sent for integration on GPU. Then, for next groups it is possible not to allocate the GPU memory and send only input data for processed elements. Also the stage of copying back the computed stiffness matrices can be eliminated in certain solution strategies, as we have mentioned already.

• 2.

  The timing results are presented for a single finite element of a given order. These results were obtained by running the kernels for a large number of elements (specified for our test runs in Tables 3 and 4) and dividing the times of performing different stages of execution by the number of elements. It should be made clear that the results depend on the total number of elements, since the times for card initialization and data transfers to and from GPU memory does not scale linearly with the number of elements.

• 3.

  The nonlinear dependence mentioned above can appear for low numbers of elements. For sufficiently high numbers, the linear dependence should emerge and the results presented in Tables 5–7 should be applicable to different problem sizes by simple scaling. The linear scaling should be observed especially for higher orders of approximation where GPU calculations, that scales linearly, take more time than card initialization and data transfers. The fact of linear scaling of the numerical integration algorithm with respect to the number of elements for sufficiently high numbers of elements has been confirmed experimentally in several studies [52, 27, 25].

• 4.

  As it was mentioned already, for CPU cores and GPU compute units we parallelize the loop over elements. If the processor resources are sufficient, than, thanks to the embarrassingly parallel character of the loop over elements in the numerical integrations algorithm, the performance of calculations should scale linearly (even in the strong sense) with the number of cores and compute units. We observed such relations for different GPUs with the NVIDIA Tesla architecture when performing computational experiments reported in [53, 34].

• 5.

  It should be noted, that both types of scaling, with respect to the number of elements and with respect to the number of cores/compute units, are much more complex when assembly process is taken into account. In such a case also the dependence of performance results on a particular problem solved (due to the different properties of the global stiffness matrix) can be observed [28, 23]. In the case of numerical integration alone, the particular problem solved is irrelevant to the performance obtained.

• 6.

  Execution times were always measured on the host side. We introduced synchronisation procedures for OpenCL events in the host code and measured elapsed times for each stage of calculations. The reason for this was to explicitly indicate the overhead associated with the subsequent stages. In recent modifications to GPU programming models it is possible to overlap different stages of execution for different kernel invocations and reduce the overheads [48].

• 7.

  Measuring elapsed time means adopting a user perspective. We do it for particular stages of calculations and also for the total execution time of kernels. The times measured on the host side were always longer than the times reported by execution environment profilers. Hence the performance that we report is the performance observed by the user, which is lower than the actual performance of the hardware.

• 8.

  Times for GPU initialization were relatively constant for all orders of approximation since they were mainly associated with creation of OpenCL memory objects in GPU memory and the total size of memory objects for each order was similar, due to the assumption that by changing the number of elements per kernel the whole available GPU memory is utilized. However, for each order of approximation the times per element are different, which is implied by the changing size of the stiffness matrix for a single element.

• 9.

  For computing performance in GFlops, only the time for kernel execution (measured on the host side) was taken into account.

• 10.

  To obtain the total execution time, we measured, in an external procedure responsible in our code for global matrix assembly, the time spent in the host code and the integration kernels.

The main results from Tables 5, 6 and 7 are illustrated in Figures 2 and 3

There are several observations that can be made based on data in Tables 5, 6 and 7:

Observation 1.

The algorithm of numerical integration for higher order finite elements has the potential for high performance execution on massively multi-core architectures. The performance obtained, above 200 GFlops, and for some orders of approximation and some variants of the algorithm above 350 GFlops for GTX580 and above 60 GFlops with several results above 150 GFlops for HD5870, can be considered high, especially for NVIDIA where it is in the range of 12%-23.5% of the peak single precision floating point performance.

Observation 2.

The performance of GTX580 GPU is much higher than that of HD5870, despite the fact that Radeon HD5870 theoretical peak performance is almost twice as big as that of GeForce GTX580. This can be explained to certain degree by higher overhead due to driver software for Radeon [27] and the fact that basic execution units of HD5870 are vector units, the fact that we do not exploit in our kernels.

Observation 3.

For GPUs, the times for global memory initialization and data transfers from the host memory are comparable to the time of kernel execution (for lower values of $p$ being even several times larger). This, once again, here for finite element numerical integration algorithm, confirms the fact that slow PCIe connection is an obstacle on a way to get higher performance for general purpose codes executed on GPUs. For the algorithm of numerical integration, with increasing order of approximation the ratio of operations to memory transfers grows and hence the influence of low PCIe performance diminishes.

Observation 4.

The times of kernel execution, i.e. the times for numerical integration calculations are always much lower for GPUs than for the reference CPU, however, the total execution times, including memory transfers are better for GPUs only in the case of higher orders of approximation.

Observation 5.

There is a large difference between the performance of the GPU as perceived by an external user and as exhibited during actual calculations. One reason for this, data transfer times, was already mentioned. The second is the fact, that, in order to ensure proper mapping of calculations onto the hardware, the GPU had to perform more operations than the CPU.

Observation 6.

There is no single kernel being the best for all situations (different orders of approximation, different hardware), although the REG_NOJAC variant is the fastest in most of cases. Different kernels exhibit similar performance, with the best results obtained for medium orders of approximation.

Observation 7.

The option of sending to GPUs, not only reference element shape functions, but also Jacobian terms for all real processed elements and all integration points (NOJAC versions) turned out to be more efficient in most cases than the option of calculating the values on GPUs (JAC variants). Despite the fact that NOJAC versions require much longer times for transferring data to GPU memory than JAC variants, this time remains much shorter than GPU initialization time and time for transfer of output data.

Observation 8.

SHM variants are usually slower than REG versions, however we believe it is worth having them both, since for different GPU architectures this situation may change.