# Parallel Implementation of FEM Solver for Shared Memory Using OpenMP

Authors: Rabia Majeed; Ramisha Farrukh; Omer Riaz; Sikandar Ali; Ali Samad; Mukhtaj Khan
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

FEM solver is designed for shared memory architectures using an object-oriented paradigm. The FEM solver is divided into five major sections; these are a collection of mesh data, creating an object for each finite element, application of boundary conditions, assembly of a linear system of equations, and lastly a solver to calculate the solution for the system of linear equations. This work has focused on the study of multiple parallel implementations for each stage of the FEM solver. All the stages are implemented in C++ while each stage parallel regions are implemented using OpenMP. The data structures used in this work are also C++ objects implemented using STL containers, like sparse matrix used is constructed as two-layer nested binary trees where each internal tree represents a sparse row in a global matrix. These structures are used in every stage and have never been converted into more conventional matrix structures like CSR, etc. Finally, this work has shown the contribution of time for each stage and the amount of speedup achieved for them. Using 16 threads, the solver has attained a speedup of 5.21 for mesh having 2099201 unknowns and Poisson’s equation. Two linear solvers Jacobi and CG are used which were consuming 95% of total solver time. The novelty of our work lies in the fact that this work focused on the implementation of Finite Element Application (FEA) using object-oriented programming for a shared memory architecture using OpenMP without considering conventional techniques of application speedup and specialized data structures like CSR, etc.


---

## Body

## 1. Introduction

Since 1999, when Intel has launched its last single-core Pentium IV processor, software development has changed drastically. Till 1999, the performance of an application was dependent on upcoming faster memory, faster processor, and memory access patterns. Due to hardware limitations, primarily heat dissipation issue Pentium IV processors were discontinued and multicore processors were introduced in the year 2006 initially as dual-core processors. The multicore processors have opened new paradigms for developers as not these have to implement their application target to utilize more than one low power core. Software engineers have to distribute their computations among available cores using either threads or processes.Finite Element Method (FEM) is a numerical method to calculate an approximate solution of ordinary differential equations and partial differential equations when provided with initial and boundary conditions. For the last decades, due to its solid theoretical foundation, it has become the most commonly used method for approximation of differential equations. FEM applications are mostly divided into four steps. These are discretization of a domain, construction of the system of equations, application of constraints, and finally approximate solution computation for the system of equations. The approximation capability of a system depends on how fine a domain is approximated. In this work, the main focus is to implement FEM applications for shared memory architectures using the OOP paradigm and Open MP library. Open MP is considered a de facto standard for multithread applications. The application developed during this research has been divided into the following four sections.


### 1.1. Read Mesh File

A mesh is an approximation of the target domain, for which calculations are performed. It is an approximate modulation of a physical object which can be anything from a high rise building to human cells. The meshes used in this work are of square shape Ω2D = (0, 1)2 and were developed using a third party software called GMSH [1]. A mesh is composed of nodes and elements where all the nodes, which have either x or y coordinates, are 0, laying on Dirichlet boundaries ΓD = 0. The provided mesh is dependent on the FE solver user and has to be read every time an FE solver is at the initial stage. For more accuracy, usually finer meshes are used which means a greater number of nodes and elements. The mesh represents a digital approximation of some real-time structure which varies from domain to domain such as a building, oceans, human bone, or cell structure. For this reason, mesh density varies due to the problem as well as error. The mesh has to change very frequently.


### 1.2. Construction of Mesh Elements

In the second stage of solver, C++ objects are created for every node as well as an element in mesh. The node objects store coordinates while the element contains store node ids, equation constants, force applied, material properties, etc. There is an associative relation between nodes and elements which enabled to sharing of node objects among element objects.


### 1.3. Construction of Linear System of Equations

Each element in Ω acts as subdomain Ωk where problem equations such as Poisson, heat, and convection diffusion are approximated. The error in solution approximation is directly proportional to the size of Ωk. The reduction in Ωk size results in a larger system of equations as more nodes are needed to construct smaller Ωk. It is always desirable to reduce the element size so that the error comes under the desired level. Global stiffness matrix A is constructed from each element contribution Ak as well as right-hand side vector b from bk.


### 1.4. Application of Dirichlet Constraints

The physical problem definition also includes different boundary conditions like Dirichlet and Neumann, etc. In this work, Dirichlet constraints are applied by FE solver application while Neumann constraints are implemented at Ωk level.


### 1.5. Solve the Linear System of Equations Using Iterative Solvers

In the last stage of this work, two iterative linear algebra solvers Jacobi and Conjugate Gradient are implemented to solve the linear system of equations.


## 1.1. Read Mesh File

A mesh is an approximation of the target domain, for which calculations are performed. It is an approximate modulation of a physical object which can be anything from a high rise building to human cells. The meshes used in this work are of square shape Ω2D = (0, 1)2 and were developed using a third party software called GMSH [1]. A mesh is composed of nodes and elements where all the nodes, which have either x or y coordinates, are 0, laying on Dirichlet boundaries ΓD = 0. The provided mesh is dependent on the FE solver user and has to be read every time an FE solver is at the initial stage. For more accuracy, usually finer meshes are used which means a greater number of nodes and elements. The mesh represents a digital approximation of some real-time structure which varies from domain to domain such as a building, oceans, human bone, or cell structure. For this reason, mesh density varies due to the problem as well as error. The mesh has to change very frequently.


## 1.2. Construction of Mesh Elements

In the second stage of solver, C++ objects are created for every node as well as an element in mesh. The node objects store coordinates while the element contains store node ids, equation constants, force applied, material properties, etc. There is an associative relation between nodes and elements which enabled to sharing of node objects among element objects.


## 1.3. Construction of Linear System of Equations

Each element in Ω acts as subdomain Ωk where problem equations such as Poisson, heat, and convection diffusion are approximated. The error in solution approximation is directly proportional to the size of Ωk. The reduction in Ωk size results in a larger system of equations as more nodes are needed to construct smaller Ωk. It is always desirable to reduce the element size so that the error comes under the desired level. Global stiffness matrix A is constructed from each element contribution Ak as well as right-hand side vector b from bk.


## 1.4. Application of Dirichlet Constraints

The physical problem definition also includes different boundary conditions like Dirichlet and Neumann, etc. In this work, Dirichlet constraints are applied by FE solver application while Neumann constraints are implemented at Ωk level.


## 1.5. Solve the Linear System of Equations Using Iterative Solvers

In the last stage of this work, two iterative linear algebra solvers Jacobi and Conjugate Gradient are implemented to solve the linear system of equations.


## 2. Background

### 2.1. Object-Oriented Finite Element Method

Finite element methods were invented to solve complex structural and elasticity problems in engineering and are frequently since the 1960s [2]. Initially, FEM solvers were implemented for supercomputers using functional programming paradigms which were cumbersome to maintain and reuse. In the 1980s, the object-oriented paradigm was the latest trend in programming practice. One of the earliest attempts to design object-oriented FEM was presented by Forde [3] in which classes were introduced for elements, materials, nodes, and boundary conditions. This work was improved by Lu et al. [4], by adding a solver class to isolate the solution algorithms from the finite element formulations. Mackie in 1992 [5] extended objects for dynamic problems to use complex geometry elements. In 1999, Mackie [6] introduced term data modeling which has enabled compatibility with other software types like CAD packages, etc. So, a class must represent data structure for basic objects; for example, degrees of freedom belong to nodes. Classes must provide access to the necessary data and services.


### 2.2. Parallel Finite Element Methods

With the advancements of finite element methods and their acceptability in the commercial sector, larger problems were designed to achieve better approximation [7]. So, the feasible means to solve such computationally intensive problems was to use parallel computers effectively [8]. Many researchers presented distributed memory designs for finite element software. One of the earliest works found was presented by K. H. Law in 1986 [9]. He proposed a solution for the MIMD machine to directly compute element distortion, without assembling the global system of equations. He purposed that, for parallel computing, iterative solution methods are suitable. Farhat and Roux in 1991 [10] presented an approach for domain decomposition regarding the solution of the parallel finite element method for equilibrium equations. The domain decomposition method was widely used earlier with supercomputers. The domain of finite elements is divided into subdomains. Each of the subdomain sets is allocated to an individual processor. Their proposed approach used fewer interprocessors communication and still offered the same amount of parallelism. Every subdomain is assigned to an individual processor and assembling of the stiffness matrix as well as the forcing vector, factoring, and computing of rigid modes have been determined. A precondition Conjugate Gradient method is applied to solve this system of equations.


### 2.3. Shared Memory Finite Element Methods

In shared memory processors, global memory (RAM) is accessed by multiple processing cores simultaneously. Since all the processors have the same view of data, these are easier to program than distributed memory processors as the communication between these processing cores is much faster as the same memory is accessed by them [11]. The same bus is used to access the main memory which limits the size of the system to only tens of processors. Sharma and Wang in 1992 [12] presented finite element programs for shared memory multiprocessors. A structure named flow node was formed which enables the programmer to make a list of blocks of computational statements and attach the properties with each relevant statement. The flow nodes are used for two purposes: an independent analysis of the code and in the derivation of numeric code for the target architecture. They generated the actual Sequent Balance code with modification [13]. Nistor et al. in 2008 [14] presented another technique of local parallelization that was to be applied on some subroutines which were CPU time-consuming inside the explicit integration main loop of the program. Only internal force vector along with stable time-step computations is parallelized by using OpenMP which generates a more efficient code without using the domain decomposition method. He presented four different techniques for parallelization. Usage of dynamic load balancing improved the efficiency of code by reducing the wait time. Turcksin et al. in 2016 [15] presented a design pattern for the multicore-enabled finite element method. In his opinion, modern codes of FEM assembling a system of linear equations are more expensive than solving them. In this regard, the authors have focused on the tasks specific to finite element toolkits. Their pattern named workstream consists of a group of elements on which local contributions are computed followed by a reduction operation. They presented three implementation schemes first of which is a parallel pipeline. Their pipeline has four stages: first is to create an item, second is to compute local contribution which can run in parallel, the third stage is the reduction, and the last stage simply marks the scratch and local contribution objects as useable again. They have not worked on only one element at a time but instead passed whole batches of elements through the pipeline. The second is the implementation of thread-local scratch elements. In the previous algorithm, scratch elements were taken round-robin from a global pool. Every processor has its cache, so this almost certainly leads to a large number of cache misses because a thread working on computing a local contribution rarely will get the same scratch objects it has while working on previous cells. So, they made the scratch objects thread-local. Local contribution objects are not thread-local. In the third implementation, they addressed another inefficiency which is the use of serial reduction operations which cannot scale to large numbers of processors if the reduction operations take a nonnegligible fraction of the runtime of the computation of the local contribution function. Reduction operations are sequential because they write into a global shared object. These memory write operations can conflict if these happen at the same time. However, in matrix assembly not all write conflict with others. For continuous elements, the reduction operations for cells conflict with each other; only these cells share degrees of freedom. So they colorize the set of cells in such a way that cells for which the reduction operations conflict are assigned different colors. With this implementation, there is no need for using create item stage because, for a given color, the order in which the elements are processed does not matter. They also presented some strategies for distributed system preceding implementations for shared memory systems. The first strategy they present is duplicating the entire data structure. This works as long as the output data structure is relatively small and the final reduction is cheap. The second approach is splitting data structures into chunks. This saves memory compared to the previous approach. But there are also downsides that load unbalancing. They used 4 16- core AMD Opteron 6378 processors running at 2.4 GHz, for a total of 64 cores to evaluate the performance. They achieved sizeable speedups between 7 and 44.


## 2.1. Object-Oriented Finite Element Method

Finite element methods were invented to solve complex structural and elasticity problems in engineering and are frequently since the 1960s [2]. Initially, FEM solvers were implemented for supercomputers using functional programming paradigms which were cumbersome to maintain and reuse. In the 1980s, the object-oriented paradigm was the latest trend in programming practice. One of the earliest attempts to design object-oriented FEM was presented by Forde [3] in which classes were introduced for elements, materials, nodes, and boundary conditions. This work was improved by Lu et al. [4], by adding a solver class to isolate the solution algorithms from the finite element formulations. Mackie in 1992 [5] extended objects for dynamic problems to use complex geometry elements. In 1999, Mackie [6] introduced term data modeling which has enabled compatibility with other software types like CAD packages, etc. So, a class must represent data structure for basic objects; for example, degrees of freedom belong to nodes. Classes must provide access to the necessary data and services.


## 2.2. Parallel Finite Element Methods

With the advancements of finite element methods and their acceptability in the commercial sector, larger problems were designed to achieve better approximation [7]. So, the feasible means to solve such computationally intensive problems was to use parallel computers effectively [8]. Many researchers presented distributed memory designs for finite element software. One of the earliest works found was presented by K. H. Law in 1986 [9]. He proposed a solution for the MIMD machine to directly compute element distortion, without assembling the global system of equations. He purposed that, for parallel computing, iterative solution methods are suitable. Farhat and Roux in 1991 [10] presented an approach for domain decomposition regarding the solution of the parallel finite element method for equilibrium equations. The domain decomposition method was widely used earlier with supercomputers. The domain of finite elements is divided into subdomains. Each of the subdomain sets is allocated to an individual processor. Their proposed approach used fewer interprocessors communication and still offered the same amount of parallelism. Every subdomain is assigned to an individual processor and assembling of the stiffness matrix as well as the forcing vector, factoring, and computing of rigid modes have been determined. A precondition Conjugate Gradient method is applied to solve this system of equations.


## 2.3. Shared Memory Finite Element Methods

In shared memory processors, global memory (RAM) is accessed by multiple processing cores simultaneously. Since all the processors have the same view of data, these are easier to program than distributed memory processors as the communication between these processing cores is much faster as the same memory is accessed by them [11]. The same bus is used to access the main memory which limits the size of the system to only tens of processors. Sharma and Wang in 1992 [12] presented finite element programs for shared memory multiprocessors. A structure named flow node was formed which enables the programmer to make a list of blocks of computational statements and attach the properties with each relevant statement. The flow nodes are used for two purposes: an independent analysis of the code and in the derivation of numeric code for the target architecture. They generated the actual Sequent Balance code with modification [13]. Nistor et al. in 2008 [14] presented another technique of local parallelization that was to be applied on some subroutines which were CPU time-consuming inside the explicit integration main loop of the program. Only internal force vector along with stable time-step computations is parallelized by using OpenMP which generates a more efficient code without using the domain decomposition method. He presented four different techniques for parallelization. Usage of dynamic load balancing improved the efficiency of code by reducing the wait time. Turcksin et al. in 2016 [15] presented a design pattern for the multicore-enabled finite element method. In his opinion, modern codes of FEM assembling a system of linear equations are more expensive than solving them. In this regard, the authors have focused on the tasks specific to finite element toolkits. Their pattern named workstream consists of a group of elements on which local contributions are computed followed by a reduction operation. They presented three implementation schemes first of which is a parallel pipeline. Their pipeline has four stages: first is to create an item, second is to compute local contribution which can run in parallel, the third stage is the reduction, and the last stage simply marks the scratch and local contribution objects as useable again. They have not worked on only one element at a time but instead passed whole batches of elements through the pipeline. The second is the implementation of thread-local scratch elements. In the previous algorithm, scratch elements were taken round-robin from a global pool. Every processor has its cache, so this almost certainly leads to a large number of cache misses because a thread working on computing a local contribution rarely will get the same scratch objects it has while working on previous cells. So, they made the scratch objects thread-local. Local contribution objects are not thread-local. In the third implementation, they addressed another inefficiency which is the use of serial reduction operations which cannot scale to large numbers of processors if the reduction operations take a nonnegligible fraction of the runtime of the computation of the local contribution function. Reduction operations are sequential because they write into a global shared object. These memory write operations can conflict if these happen at the same time. However, in matrix assembly not all write conflict with others. For continuous elements, the reduction operations for cells conflict with each other; only these cells share degrees of freedom. So they colorize the set of cells in such a way that cells for which the reduction operations conflict are assigned different colors. With this implementation, there is no need for using create item stage because, for a given color, the order in which the elements are processed does not matter. They also presented some strategies for distributed system preceding implementations for shared memory systems. The first strategy they present is duplicating the entire data structure. This works as long as the output data structure is relatively small and the final reduction is cheap. The second approach is splitting data structures into chunks. This saves memory compared to the previous approach. But there are also downsides that load unbalancing. They used 4 16- core AMD Opteron 6378 processors running at 2.4 GHz, for a total of 64 cores to evaluate the performance. They achieved sizeable speedups between 7 and 44.


## 3. FE Solver

Figure1 shows a class diagram of implemented Finite Element Poisson 2D solver. Poisson_2D_Solver class is the main class that has four public methods that represent each stage of finite element solutions. The generateElements method is used to read the mesh file and create objects for nodes and elements. The node class is implemented to store each node's ids and coordinates (x, y). An abstract element class is represented to define a common interface for all kinds of subdomains. The Poisson 2D element classes Edge and Triangle are inherited from element class and their objects will be accessed from the element pointer using polymorphism. Gmsh_Reader class is implemented to read nodes and elements data one by one using getNode and getElement methods, respectively.Figure 1
Class diagram of FE Poisson solver.Poisson_2D_Solver::generateEquation method constructs a system of equations using mesh element objects. In this function, all the elementsAk and bk are calculated and accumulated into A=∑ΛkAkΛkT and b=∑Λkbk, where Ak and bk are Ωk contributions to global system matrix A and load vector b while Λ is used to map local node’s ids to global node ids.Each element is composed of a unique subset of nodes, and this information is private to every element. There are many compressed sparse matrix formats like CSR (Compressed Sparse Row), COO (Coordinate Format), BSR (Block Sparse Row), etc. for efficient sparse matrix-vector multiplication. These formats need to know in advance the number of nonzero entries in each row and their column indexes. The mapping among nodes is not known in this interface. A dynamic sparse matrix data structure is required, which should have constant time for insertion, searching, and deletion. Sparse_Matrix class is implemented to storeA, which has a complexity of logn. The Vector class is implemented to store b with a complexity of 1.Poisson_2D_Solver::applyDirichlet method applies Dirichlet constraints onA and b. The list of Dirichlet nodes is collected during the elements’ generation stage. The user has to provide information on which boundary of Ω is Dirichlet and the amount of force is being applied on these. In this work, all the nodes which have either x or y coordinate or both are zero lays at Dirichlet boundary and zero force is being applied on Dirichlet boundaries.Poisson_2D_Solver::getSolution method calculates the solution of the linear system of equations in which two iterative linear solvers Jacobi and Conjugate Gradient are implemented.


## 4. Parallel FE Solver

All stages mentioned above are implemented in parallel for shared memory architecture using OpenMP. The selection of OpenMP was because it is a directive-based language that is very convenient to adopt and secondly software market has opted for it as a de facto standard. It is implemented by all major C++ compilers like Microsoft, Intel, GNU, etc. OpenMP hides the details about creating, management, and deletion of threads. A developer only indicates where to create a parallel region and which elements have to pass a shared or private variable. A developer also has to identify critical sections and methods to manage these critical sections.


### 4.1. Read Mesh

There were multiple algorithms considered during the parallel generateElement method.


#### 4.1.1. Single File Handler Shared among Threads

The first approach considered for parallel mesh read is to share a file handler among multiple threads. Initially, the file is read in a critical section where it is observed that mesh read time increases with a rise in the number of threads due to the fact all threads except one are waiting for the critical section. Secondly, a shared file handler is used in each thread without a critical section. It resulted in performance gain but it was observed that some nodes and elements are read multiple times and their objects are created more than one time.


#### 4.1.2. Each Thread Has Its Own File Handler

The second approach considered provides a private file handler to each thread. All threads initially calculate how many nodes and elements each thread will read and what should be their starting ids and end ids. Figure2 shows how nodes and elements are mapped to threads using block distribution and cyclic distribution, respectively. Block distribution allocates consecutive rows to each thread as shown in Figure 2(a) while cyclic distribution allocates every ith row to each thread where there is the number of threads as shown in Figure 2(b). Initial experiments showed that block distribution has outperformed cyclic distribution as can be seen in Tables 1 and 2. In both distributions, an almost equal number of nodes and elements are read while each node and element are read in a single parallel region only. For 16 threads, cyclic distribution took 9.84 secs while block distribution took 1.24 seconds.Figure 2
Block distribution and cyclic distribution. (a) Block distribution. (b) Cycle distribution.
(a)(b)Table 1
Speedup and efficiency of mesh read module using cyclic distribution (525313 and 1050628).
Threads12345678Speedup10.861.131.151.211.311.091.03Efficiency10.430.380.290.240.220.160.13Threads910111213141516Speedup0.990.900.870.840.810.800.790.76Efficiency0.110.090.080.070.060.060.050.05Table 2
Speedup and efficiency of mesh read module using block distribution (525313 nodes and 1050628 elements).
Threads12345678Speedup1.0001.8302.7383.5254.1874.8895.5206.089Efficiency100.00091.48591.27188.13783.73781.48378.85776.109Threads910111213141516Speedup6.6517.1467.5848.1408.3108.8039.1009.107Efficiency73.89571.46568.94367.83763.92162.88160.66756.917


### 4.2. Construction of Equations

Constructing a system of equations from a mesh requires contribution from each element of mesh. The system of the equation is composed of stiffness matrixA, solution vector u, and load vector b. For each element k, local stiffness matrix Ak and load vector bk are computed and added to A and b, respectively. In (1) and (2), Λk is the connectivity matrix used to introduce the local element k system of equations to the global system of equations.(1)A=∑ΛkAkΛkT,(2)b=∑Λkbk.In the parallel version, multiple approaches are adopted, which are discussed below with their advantages and disadvantages. The application was designed in a way to use less memory footprint as a well efficient method of calculation. Element class interface which is inherited by all element classes (as can be seen in Figure1) provides a reference of global structures to each element while allowing each element to directly store local data into global data structures. In the parallel implementation, multiple elements add their contribution to global elements simultaneously which makes this piece of code a source of error. This critical section was addressed using the following methods.


#### 4.2.1. Thread-Safe Data Structures

In this approach, a thread-safe data structure forA and b is considered. In the first structure of A, shown in Figure 3(a), a single lock is added to sparse_matrix class and each time any of the threads either performs data insertion, modification. or deletion, lock must be acquired. The acquisition and release of lock are added to data structure insertion and deletion methods. This approach resulted in correct matrix data but the time taken was similar to the sequential time as shown in Tables 3 and 4. In both cases, maximum speedup was attained for three threads and reduces as threads increased. Only one thread has access to the critical section while other threads are waiting.Figure 3
Locking mechanisms adopted for the sparse matrix data structure. (a) Single lock for a matrix. (b) Separate lock for each row.
(a)(b)Table 3
Construction of equation using global locking using mesh (525313 nodes and 1050628 elements).
Threads12345678Time9.10044.9133.7544.7784.9975.1215.8936.735Speedup11.8522.4241.9051.8221.7771.5441.351Threads910111213141516Time6.3717.0626.9148.2088.4567.32710.5128.186Speedup1.4281.2891.3161.1091.0761.2420.8661.112Table 4
Construction of equation using global locking using mesh (2099201 nodes and 4198404 elements).
Threads12345678Time52.2725.81618.12120.04222.94423.60321.97422.677Speedup12.0252.8862.6082.2782.2152.3792.305Threads910111213141516Time21.30126.98330.51140.22041.60634.91247.10745.224Speedup2.4541.9371.7131.2991.2561.4981.1091.156It is observed that each element is constructed using a unique subset of nodes which means that not all the elements are updating the same rows at any instance. The sparse_matrix is composed of two-level nested binary trees using standard template library binary tree implementation using std:map provided in Standard Template Library (STL). The sparse_matrix class is modified so that a lock is created for each row as shown in Figure3(b). So, only those threads which are modifying the same rows have to wait. Tables 5 and 6 show the timing of the row-level locking mechanism. These have shown a continuous rise in speedup with an increase in threads as in most cases threads will be writing to different rows and very few threads are waiting at any time.Table 5
Construction of equation using row locking using mash (525313 nodes and 1050628 elements).
Threads12345678Time12.5975.8953.6413.1462.4772.2151.8031.750Speedup1.0002.1373.4604.0045.0865.6886.9897.198Threads910111213141516Time1.6201.3121.3361.1841.0771.0670.9790.898Speedup7.7749.6039.43110.64011.69911.81112.87214.025Table 6
Construction of equation using row locking using mash (2099201 nodes and 4198404 elements).
Threads12345678Time63.6728.4419.1114.7512.2611.039.528.82Speedup1.002.243.334.325.205.786.697.22Threads910111213141516Time7.787.106.376.235.715.765.145.02Speedup8.188.9710.0010.2311.1511.0612.3912.67


#### 4.2.2. Use OpenMP Critical Section to Gather Equation Data

The second approach was to implement a critical section at the application in which a non-thread-safe data structure is used forA and b. The parallel region was developed using #pragma omp for where each element updates the global data structure in a critical section. The timing of this implementation is shown in Table 7. This method has shown no speedup at all as most of the time, except for one thread, all are in the waiting stage. Even though the number of tasks per thread reduces with a rise in the number of threads but wait time also increases.Table 7
Construction of equation using OpenMP critical using mash (2099201 nodes and 4198404 elements).
Threads12345678Time61.9744.9253.6048.4246.5250.4548.4750.02Speedup1.001.381.161.281.331.231.281.24Threads910111213141516Time48.6254.3350.3450.9051.1556.3552.4356.74Speedup1.281.141.231.221.211.101.181.09


#### 4.2.3. Use Threads’ Local Data Structures and Gather Using Butterfly Communication

This approach aims to void critical sections by using the local sparse matrix in each thread. As shown in Figure4, this algorithm is composed of two sections where in earlier stage each thread collects data from each of its designated elements into local data structures while the latter stage is to iteratively add these threads’ local data structures At and bt to A and b using butterfly communications. The butterfly communication is an iterative process where the number of iterations is equal to log2(total threads). Figure 4 has only log2(8) = 3 synchronization shown by an orange horizontal line. Tables 8 and 9 show timing where each table has three rows. The first row indicates the construction of At by summing all the allocated elements’ stiffness matrices. The second row indicates the time for the construction of a global matrix A from At while the last row indicates speedup achieved by multiple threads. Each thread has to sequentially construct its At. The population of A requires data from At using butterfly communication. More matrices are added as the number of threads increases and more stages are added to butterfly communication. During experiments, elements are read in the same order as given in the file and a nonoverlapping subset of consecutive elements was allocated to each thread. Finally, the increase in the number of threads reduces At and increases time for A. Overall time remains constant and speedup always remained under 2%.Figure 4
Construction of equation using butterfly communications.Table 8
Construction of equation using butterfly communication using mash (525313 nodes and 1050628 elements).
Threads12345678At=∑Ak22.43110.6256.8995.0493.9703.2572.8182.437A=∑At05.5577.2448.2009.0199.0549.4999.541Speedup1.0001.3771.4951.6611.6981.8071.8111.843Threads910111213141516At=∑Ak2.1651.9521.7161.5981.4761.3841.2831.186A=∑At10.39510.11310.33910.58110.47310.66810.52510.577Speedup1.7741.8491.8211.8071.8661.8531.8791.849Table 9
Construction of equation using butterfly communication using mash (2099201 nodes and 4198404 elements).
Threads12345678At=∑Ak99.435747.37130.87122.81317.41414.91712.42910.996A=∑At024.49633.24937.38139.70741.88742.69244.496Speedup1.0001.3841.5091.6511.7151.7451.7901.789Threads910111213141516At=∑Ak9.8548.8617.9627.2966.6516.0825.8035.486A=∑At45.82945.24944.60845.31745.63245.67446.64947.515Speedup1.7801.8201.8821.8861.8901.9081.8881.868


### 4.3. Application of Constraints

In this work, only the Dirichlet constraint is considered. The parallel algorithm for Dirichlet constraints has two parts as shown in Figure5. The first part is to remove all the Dirichlet rows from A while setting its diagonal entry to 1. In this experiment, a two-level nested std:map is used to store a sparse matrix. The std:map is not a thread-safe data structure. The removal of rows means the removal of internal std:map object from outer std:map object where each edge represents a row. The second part is more time-consuming as for every non-Dirichlet row, Dirichlet columns entries are searched and removed. Figure 6 shows a 7 × 7 stiffness matrix before and after applying Dirichlet constraints where nodes 4 and 6 lay on the Dirichlet boundary. Removal of Dirichlet rows was performed sequentially because of a non-thread-safe data structure. Non-Diriechlet rows were distributed among the threads where each thread will remove Dirichlet columns from each row. Table 10 shows timing and speedup for the Dirichlet algorithm. The second part of the algorithm is embarrassingly parallel and has achieved an almost linear speedup of 14.65%for 16 threads.Figure 5
Parallel Dirichlet implementation.Figure 6
Dirichlet constraints. (a) Global stiffness matrix. (b) Global stiffness matrix after applying Dirichlet constraints.
(a)(b)Table 10
Dirichlet application using multiple threads using a problem having 2099201 unknown.
Threads12345678Time1671825551423343281241218Speedup1.002.033.033.954.875.966.927.68Threads910111213141516Time194173158149138128119114Speedup8.619.6810.5811.2512.1613.1114.0614.65


### 4.4. Linear Algebra Solvers

There were two most commonly used iterative solvers Jacobi and Conjugate Gradient considered in this work. These are implemented in parallel using OpenMP. Their algorithms and implementation details are as follows. All the critical sections are colored blue while synchronization points are colored red.(1)
Jacobi solver iteration can be segregated into the following parts.(i)
Calculation of residualr=b−Ax(ii)
New approximation ofx+=D−1r(iii)
Calculation of error thresholde=r/x(2)
Conjugate Gradient SolverAlgorithm 1: Jacobi solver.
Input:A, b⟶,x, N, x_norm, r_normOutput:x(1)
fork (A, b⟶,x, N, x_norm, r_norm)(2)
Th_id, Tt_id (Thread id and Total Thread by OS)(3)
St_RowId ← Ed_RowId ← 0 {Start and end ids}(4)
IfTh_id ＜ N mod Tt_ththen(5)
St_RowID ← Th_id∗N/Tt/th+Th_id(6)
Ed_RowID ← St_RowID+N/Tt_th+1(7)
else(8)
St_RowID ← Th_id ∗ N/Tt_th+Tt_th(9)
Ed_RowID ← St_RowID+N/Tt_th(10)
end_if(11)
r {Local Residual vector for allocated nodes}(12)
error ← 1(13)
Whileerror > Thresholddo(14)
r_tmp ← x_tmp ← 0 {Local Norms for thread}(15)
foreachi ← St_RowID, Ed,RowIDdo(16)
repeat(17)
r[i-St_RowID]←b[i]-A[i,j]∗x[j](18)
untilFor each entry j for row i of A(19)
r_tmp ← r_tmp + r[i-St_RowID]2(20)
end for(21)
x_norm ← r_norm ← 0 [Global Norms](22)
barrier(23)
foreachi ← St_RowID, Ed_RowIDdo(24)
x[i] ← x[i] + r[i-St_RowID]/A[i,i](25)
x_tmp ← x_tmp+x[i]2(26)
end for(27)
lock(28)
r_norm ←r_norm+r_tmp(29)
x_norm←x_norm+x_tmp(30)
unlock(31)
barrier(32)
error ←r_norm/x_norm(33)
end while(34)
joinFigure 7
Parallel Jacobi solver stages. (a) Parallel residual calculation. (b) Error calculation.The second implementation for Jacobi is to create a single parallel region that will cover all three parts of Jacobi iteration as shown in Algorithm1. There is only one parallel region in this application starting at line 1 and collapsing at line 35. All three stages are implemented in single parallel regions, and different stages have to be segregated using explicit barriers at line 22. The first explicit barrier is necessary as in the first region x is used to calculate r and in the second region x is updated using r. The second explicit barrier at line 31 is required to the calculation of error which requires r/x. Before the second barrier, a critical section, from lines 27 to 30, is required to calculate r and x. Table 11 shows timing data for Algorithm 1 calculated using a mesh of 131585 nodes and 263168 elements. Jacobi has a slow convergence rate that is taking a large number of iterations to converge and it is the most time-consuming part of FE solver. For this problem, it took about 165009 iterations. The parallel task has enabled the reduction of the Jacobi solver's sequential time of 3780 sec to parallel time of 605 sec using 16 threads and achieving 6.246 times speedup. During execution, the stopping criteria adopted was either the error is less than 10–8 or total iterations are greater than 106.Table 11
Jacobi solver parallel timing and speedup for system of equations having 131585 unknowns.
Threads12345678Jacobi3780258320231458145511781129982Speedup1.001.461.872.592.603.213.353.85Threads910111213141516Jacobi899814754733681642619605Speedup4.204.645.015.165.555.896.116.25Algorithm2 shows a parallel algorithm of Conjugate Gradient. While designing this algorithm, the main aim was to keep the global synchronization minimum. In this algorithm, global synchronization is represented by barrier function at lines 17, 25, and 34. All the codes displayed in the black color are each thread individual code that has to perform without communicating with other threads. All the codes in blue are a critical section that has to be executed by all threads but one at a time. In this section, global variables rtr1, rtr2, and ptAp are computed. In each iteration, there are two synchronization points required to update shared variables such as ptAp and rtr2, respectively.Algorithm 2: Parallel Conjugate Gradient algorithm.
Input: (A, b⟶,x, N, rtt1, rtr2, ptAp)Output:x(1)
fork (A, b⟶,x, N, rtt1, rtr2, ptAp) {shared variables}(2)
Th_id, Tt_id (Thread id and Total Thread by OS)(3)
St_RowId ← Ed_RowId ← 0 {Start and end ids}(4)
IfTh_id ＜ N % Tt_ththen(5)
St_RowID ← Th_id∗N/Tt+Th_id(6)
Ed_RowID ← St_RowID+N/Tt_th+1(7)
else(8)
St_RowID ← Th_id ∗ N/Tt_th+Tt_th(9)
Ed_RowID ← St_RowID+N/Tt_th(10)
end_if(11)
rtr1 ← rtr2 ← t_rtr ← 0(12)
foreachi ← St_RowID.Ed_RowIDdo(13)
p[i] ← r[i] ← b[i]-Ai:,x(14)
t_rtr ← t_rtr + r[i] ∗ r[i](15)
end for(16)
lockrtr2 ← rtr2 + t_rtr unlock(17)
barrier(18)
whilertr2>Threshold do(19)
rtr1 ← rtr2(20)
rtr2 ← t_ptAp ← 0(21)
foreachi ← St_RowID, Ed_RowIDdo(22)
t_p'Ap ← t_ptAp+p[i].Ai:,p(23)
end for(24)
lockptAp ← ptAp+t_ptAp unlock(25)
barrier(26)
t_rtr ← 0(27)
foreachi ← St_RowID, Ed_RowIDdo(28)
x[i] ← x[i]+(rtr1/ptAp).p[i](29)
r[i] ← r[i]+(rtr1/ptAp).Ai:,p(30)
t_rtr ← t_rtr+r[i] ∗ r [i](31)
end for(32)
ptAp ← 0(33)
lockrtr2 ← rtr2+t_rtr unlock(34)
barrier(35)
foreachi ← St_RowID, Ed_RowIDdo(36)
p[i] ← r[i]+(rtr2/rtr1).p[i](37)
end for(38)
end while(39)
joinTables12 and 13 show the timing and speedup data of Algorithm 2 for two mashes. The system of 131585 unknown took 1019 iterations while the system of 2099281 unknown took 3987 iterations. Algorithm 2 has three critical sections (all colored in blue) where the first critical section at line 16 has to be executed once as it is before while loop while the other two sections at lines 24 and 33 are in a while loop and have to be executed in every iteration. These sections are to update global variables (rtr and ptAp) with thread-local calculations. These variables cannot be updated in the same critical section as ptAp to use for the calculation of rtr. Therefore, every time a thread after going through a critical section has to wait for other threads to execute it before further processing, this is imposed by global communication. For a fixed size problem, as the number of threads in the parallel region increases, the time required for thread-independent tasks (colored in black grey) will reduce while the critical section and communication timing will increase.Table 12
Conjugate Gradient solver parallel timing and speedup for system of equations having 131585 unknowns.
Threads12345678CG23.4217.2612.6410.108.437.296.906.10Speedup1.001.361.852.322.783.213.393.85Threads910111213141516CG5.465.124.694.394.233.853.653.62Speedup4.294.575.005.345.546.096.426.47Table 13
Conjugate Gradient solver parallel timing and speedup for system of equations having 2099201 unknowns.
Threads12345678CG15661214997931721687670634Speedup1.001.291.571.682.172.282.342.47Threads910111213141516CG520519.95479449399366.93366312Speedup3.013.013.273.493.934.274.285.01


## 4.1. Read Mesh

There were multiple algorithms considered during the parallel generateElement method.


### 4.1.1. Single File Handler Shared among Threads

The first approach considered for parallel mesh read is to share a file handler among multiple threads. Initially, the file is read in a critical section where it is observed that mesh read time increases with a rise in the number of threads due to the fact all threads except one are waiting for the critical section. Secondly, a shared file handler is used in each thread without a critical section. It resulted in performance gain but it was observed that some nodes and elements are read multiple times and their objects are created more than one time.


### 4.1.2. Each Thread Has Its Own File Handler

The second approach considered provides a private file handler to each thread. All threads initially calculate how many nodes and elements each thread will read and what should be their starting ids and end ids. Figure2 shows how nodes and elements are mapped to threads using block distribution and cyclic distribution, respectively. Block distribution allocates consecutive rows to each thread as shown in Figure 2(a) while cyclic distribution allocates every ith row to each thread where there is the number of threads as shown in Figure 2(b). Initial experiments showed that block distribution has outperformed cyclic distribution as can be seen in Tables 1 and 2. In both distributions, an almost equal number of nodes and elements are read while each node and element are read in a single parallel region only. For 16 threads, cyclic distribution took 9.84 secs while block distribution took 1.24 seconds.Figure 2
Block distribution and cyclic distribution. (a) Block distribution. (b) Cycle distribution.
(a)(b)Table 1
Speedup and efficiency of mesh read module using cyclic distribution (525313 and 1050628).
Threads12345678Speedup10.861.131.151.211.311.091.03Efficiency10.430.380.290.240.220.160.13Threads910111213141516Speedup0.990.900.870.840.810.800.790.76Efficiency0.110.090.080.070.060.060.050.05Table 2
Speedup and efficiency of mesh read module using block distribution (525313 nodes and 1050628 elements).
Threads12345678Speedup1.0001.8302.7383.5254.1874.8895.5206.089Efficiency100.00091.48591.27188.13783.73781.48378.85776.109Threads910111213141516Speedup6.6517.1467.5848.1408.3108.8039.1009.107Efficiency73.89571.46568.94367.83763.92162.88160.66756.917


## 4.1.1. Single File Handler Shared among Threads

The first approach considered for parallel mesh read is to share a file handler among multiple threads. Initially, the file is read in a critical section where it is observed that mesh read time increases with a rise in the number of threads due to the fact all threads except one are waiting for the critical section. Secondly, a shared file handler is used in each thread without a critical section. It resulted in performance gain but it was observed that some nodes and elements are read multiple times and their objects are created more than one time.


## 4.1.2. Each Thread Has Its Own File Handler

The second approach considered provides a private file handler to each thread. All threads initially calculate how many nodes and elements each thread will read and what should be their starting ids and end ids. Figure2 shows how nodes and elements are mapped to threads using block distribution and cyclic distribution, respectively. Block distribution allocates consecutive rows to each thread as shown in Figure 2(a) while cyclic distribution allocates every ith row to each thread where there is the number of threads as shown in Figure 2(b). Initial experiments showed that block distribution has outperformed cyclic distribution as can be seen in Tables 1 and 2. In both distributions, an almost equal number of nodes and elements are read while each node and element are read in a single parallel region only. For 16 threads, cyclic distribution took 9.84 secs while block distribution took 1.24 seconds.Figure 2
Block distribution and cyclic distribution. (a) Block distribution. (b) Cycle distribution.
(a)(b)Table 1
Speedup and efficiency of mesh read module using cyclic distribution (525313 and 1050628).
Threads12345678Speedup10.861.131.151.211.311.091.03Efficiency10.430.380.290.240.220.160.13Threads910111213141516Speedup0.990.900.870.840.810.800.790.76Efficiency0.110.090.080.070.060.060.050.05Table 2
Speedup and efficiency of mesh read module using block distribution (525313 nodes and 1050628 elements).
Threads12345678Speedup1.0001.8302.7383.5254.1874.8895.5206.089Efficiency100.00091.48591.27188.13783.73781.48378.85776.109Threads910111213141516Speedup6.6517.1467.5848.1408.3108.8039.1009.107Efficiency73.89571.46568.94367.83763.92162.88160.66756.917


## 4.2. Construction of Equations

Constructing a system of equations from a mesh requires contribution from each element of mesh. The system of the equation is composed of stiffness matrixA, solution vector u, and load vector b. For each element k, local stiffness matrix Ak and load vector bk are computed and added to A and b, respectively. In (1) and (2), Λk is the connectivity matrix used to introduce the local element k system of equations to the global system of equations.(1)A=∑ΛkAkΛkT,(2)b=∑Λkbk.In the parallel version, multiple approaches are adopted, which are discussed below with their advantages and disadvantages. The application was designed in a way to use less memory footprint as a well efficient method of calculation. Element class interface which is inherited by all element classes (as can be seen in Figure1) provides a reference of global structures to each element while allowing each element to directly store local data into global data structures. In the parallel implementation, multiple elements add their contribution to global elements simultaneously which makes this piece of code a source of error. This critical section was addressed using the following methods.


### 4.2.1. Thread-Safe Data Structures

In this approach, a thread-safe data structure forA and b is considered. In the first structure of A, shown in Figure 3(a), a single lock is added to sparse_matrix class and each time any of the threads either performs data insertion, modification. or deletion, lock must be acquired. The acquisition and release of lock are added to data structure insertion and deletion methods. This approach resulted in correct matrix data but the time taken was similar to the sequential time as shown in Tables 3 and 4. In both cases, maximum speedup was attained for three threads and reduces as threads increased. Only one thread has access to the critical section while other threads are waiting.Figure 3
Locking mechanisms adopted for the sparse matrix data structure. (a) Single lock for a matrix. (b) Separate lock for each row.
(a)(b)Table 3
Construction of equation using global locking using mesh (525313 nodes and 1050628 elements).
Threads12345678Time9.10044.9133.7544.7784.9975.1215.8936.735Speedup11.8522.4241.9051.8221.7771.5441.351Threads910111213141516Time6.3717.0626.9148.2088.4567.32710.5128.186Speedup1.4281.2891.3161.1091.0761.2420.8661.112Table 4
Construction of equation using global locking using mesh (2099201 nodes and 4198404 elements).
Threads12345678Time52.2725.81618.12120.04222.94423.60321.97422.677Speedup12.0252.8862.6082.2782.2152.3792.305Threads910111213141516Time21.30126.98330.51140.22041.60634.91247.10745.224Speedup2.4541.9371.7131.2991.2561.4981.1091.156It is observed that each element is constructed using a unique subset of nodes which means that not all the elements are updating the same rows at any instance. The sparse_matrix is composed of two-level nested binary trees using standard template library binary tree implementation using std:map provided in Standard Template Library (STL). The sparse_matrix class is modified so that a lock is created for each row as shown in Figure3(b). So, only those threads which are modifying the same rows have to wait. Tables 5 and 6 show the timing of the row-level locking mechanism. These have shown a continuous rise in speedup with an increase in threads as in most cases threads will be writing to different rows and very few threads are waiting at any time.Table 5
Construction of equation using row locking using mash (525313 nodes and 1050628 elements).
Threads12345678Time12.5975.8953.6413.1462.4772.2151.8031.750Speedup1.0002.1373.4604.0045.0865.6886.9897.198Threads910111213141516Time1.6201.3121.3361.1841.0771.0670.9790.898Speedup7.7749.6039.43110.64011.69911.81112.87214.025Table 6
Construction of equation using row locking using mash (2099201 nodes and 4198404 elements).
Threads12345678Time63.6728.4419.1114.7512.2611.039.528.82Speedup1.002.243.334.325.205.786.697.22Threads910111213141516Time7.787.106.376.235.715.765.145.02Speedup8.188.9710.0010.2311.1511.0612.3912.67


### 4.2.2. Use OpenMP Critical Section to Gather Equation Data

The second approach was to implement a critical section at the application in which a non-thread-safe data structure is used forA and b. The parallel region was developed using #pragma omp for where each element updates the global data structure in a critical section. The timing of this implementation is shown in Table 7. This method has shown no speedup at all as most of the time, except for one thread, all are in the waiting stage. Even though the number of tasks per thread reduces with a rise in the number of threads but wait time also increases.Table 7
Construction of equation using OpenMP critical using mash (2099201 nodes and 4198404 elements).
Threads12345678Time61.9744.9253.6048.4246.5250.4548.4750.02Speedup1.001.381.161.281.331.231.281.24Threads910111213141516Time48.6254.3350.3450.9051.1556.3552.4356.74Speedup1.281.141.231.221.211.101.181.09


### 4.2.3. Use Threads’ Local Data Structures and Gather Using Butterfly Communication

This approach aims to void critical sections by using the local sparse matrix in each thread. As shown in Figure4, this algorithm is composed of two sections where in earlier stage each thread collects data from each of its designated elements into local data structures while the latter stage is to iteratively add these threads’ local data structures At and bt to A and b using butterfly communications. The butterfly communication is an iterative process where the number of iterations is equal to log2(total threads). Figure 4 has only log2(8) = 3 synchronization shown by an orange horizontal line. Tables 8 and 9 show timing where each table has three rows. The first row indicates the construction of At by summing all the allocated elements’ stiffness matrices. The second row indicates the time for the construction of a global matrix A from At while the last row indicates speedup achieved by multiple threads. Each thread has to sequentially construct its At. The population of A requires data from At using butterfly communication. More matrices are added as the number of threads increases and more stages are added to butterfly communication. During experiments, elements are read in the same order as given in the file and a nonoverlapping subset of consecutive elements was allocated to each thread. Finally, the increase in the number of threads reduces At and increases time for A. Overall time remains constant and speedup always remained under 2%.Figure 4
Construction of equation using butterfly communications.Table 8
Construction of equation using butterfly communication using mash (525313 nodes and 1050628 elements).
Threads12345678At=∑Ak22.43110.6256.8995.0493.9703.2572.8182.437A=∑At05.5577.2448.2009.0199.0549.4999.541Speedup1.0001.3771.4951.6611.6981.8071.8111.843Threads910111213141516At=∑Ak2.1651.9521.7161.5981.4761.3841.2831.186A=∑At10.39510.11310.33910.58110.47310.66810.52510.577Speedup1.7741.8491.8211.8071.8661.8531.8791.849Table 9
Construction of equation using butterfly communication using mash (2099201 nodes and 4198404 elements).
Threads12345678At=∑Ak99.435747.37130.87122.81317.41414.91712.42910.996A=∑At024.49633.24937.38139.70741.88742.69244.496Speedup1.0001.3841.5091.6511.7151.7451.7901.789Threads910111213141516At=∑Ak9.8548.8617.9627.2966.6516.0825.8035.486A=∑At45.82945.24944.60845.31745.63245.67446.64947.515Speedup1.7801.8201.8821.8861.8901.9081.8881.868


## 4.2.1. Thread-Safe Data Structures

In this approach, a thread-safe data structure forA and b is considered. In the first structure of A, shown in Figure 3(a), a single lock is added to sparse_matrix class and each time any of the threads either performs data insertion, modification. or deletion, lock must be acquired. The acquisition and release of lock are added to data structure insertion and deletion methods. This approach resulted in correct matrix data but the time taken was similar to the sequential time as shown in Tables 3 and 4. In both cases, maximum speedup was attained for three threads and reduces as threads increased. Only one thread has access to the critical section while other threads are waiting.Figure 3
Locking mechanisms adopted for the sparse matrix data structure. (a) Single lock for a matrix. (b) Separate lock for each row.
(a)(b)Table 3
Construction of equation using global locking using mesh (525313 nodes and 1050628 elements).
Threads12345678Time9.10044.9133.7544.7784.9975.1215.8936.735Speedup11.8522.4241.9051.8221.7771.5441.351Threads910111213141516Time6.3717.0626.9148.2088.4567.32710.5128.186Speedup1.4281.2891.3161.1091.0761.2420.8661.112Table 4
Construction of equation using global locking using mesh (2099201 nodes and 4198404 elements).
Threads12345678Time52.2725.81618.12120.04222.94423.60321.97422.677Speedup12.0252.8862.6082.2782.2152.3792.305Threads910111213141516Time21.30126.98330.51140.22041.60634.91247.10745.224Speedup2.4541.9371.7131.2991.2561.4981.1091.156It is observed that each element is constructed using a unique subset of nodes which means that not all the elements are updating the same rows at any instance. The sparse_matrix is composed of two-level nested binary trees using standard template library binary tree implementation using std:map provided in Standard Template Library (STL). The sparse_matrix class is modified so that a lock is created for each row as shown in Figure3(b). So, only those threads which are modifying the same rows have to wait. Tables 5 and 6 show the timing of the row-level locking mechanism. These have shown a continuous rise in speedup with an increase in threads as in most cases threads will be writing to different rows and very few threads are waiting at any time.Table 5
Construction of equation using row locking using mash (525313 nodes and 1050628 elements).
Threads12345678Time12.5975.8953.6413.1462.4772.2151.8031.750Speedup1.0002.1373.4604.0045.0865.6886.9897.198Threads910111213141516Time1.6201.3121.3361.1841.0771.0670.9790.898Speedup7.7749.6039.43110.64011.69911.81112.87214.025Table 6
Construction of equation using row locking using mash (2099201 nodes and 4198404 elements).
Threads12345678Time63.6728.4419.1114.7512.2611.039.528.82Speedup1.002.243.334.325.205.786.697.22Threads910111213141516Time7.787.106.376.235.715.765.145.02Speedup8.188.9710.0010.2311.1511.0612.3912.67


## 4.2.2. Use OpenMP Critical Section to Gather Equation Data

The second approach was to implement a critical section at the application in which a non-thread-safe data structure is used forA and b. The parallel region was developed using #pragma omp for where each element updates the global data structure in a critical section. The timing of this implementation is shown in Table 7. This method has shown no speedup at all as most of the time, except for one thread, all are in the waiting stage. Even though the number of tasks per thread reduces with a rise in the number of threads but wait time also increases.Table 7
Construction of equation using OpenMP critical using mash (2099201 nodes and 4198404 elements).
Threads12345678Time61.9744.9253.6048.4246.5250.4548.4750.02Speedup1.001.381.161.281.331.231.281.24Threads910111213141516Time48.6254.3350.3450.9051.1556.3552.4356.74Speedup1.281.141.231.221.211.101.181.09


## 4.2.3. Use Threads’ Local Data Structures and Gather Using Butterfly Communication

This approach aims to void critical sections by using the local sparse matrix in each thread. As shown in Figure4, this algorithm is composed of two sections where in earlier stage each thread collects data from each of its designated elements into local data structures while the latter stage is to iteratively add these threads’ local data structures At and bt to A and b using butterfly communications. The butterfly communication is an iterative process where the number of iterations is equal to log2(total threads). Figure 4 has only log2(8) = 3 synchronization shown by an orange horizontal line. Tables 8 and 9 show timing where each table has three rows. The first row indicates the construction of At by summing all the allocated elements’ stiffness matrices. The second row indicates the time for the construction of a global matrix A from At while the last row indicates speedup achieved by multiple threads. Each thread has to sequentially construct its At. The population of A requires data from At using butterfly communication. More matrices are added as the number of threads increases and more stages are added to butterfly communication. During experiments, elements are read in the same order as given in the file and a nonoverlapping subset of consecutive elements was allocated to each thread. Finally, the increase in the number of threads reduces At and increases time for A. Overall time remains constant and speedup always remained under 2%.Figure 4
Construction of equation using butterfly communications.Table 8
Construction of equation using butterfly communication using mash (525313 nodes and 1050628 elements).
Threads12345678At=∑Ak22.43110.6256.8995.0493.9703.2572.8182.437A=∑At05.5577.2448.2009.0199.0549.4999.541Speedup1.0001.3771.4951.6611.6981.8071.8111.843Threads910111213141516At=∑Ak2.1651.9521.7161.5981.4761.3841.2831.186A=∑At10.39510.11310.33910.58110.47310.66810.52510.577Speedup1.7741.8491.8211.8071.8661.8531.8791.849Table 9
Construction of equation using butterfly communication using mash (2099201 nodes and 4198404 elements).
Threads12345678At=∑Ak99.435747.37130.87122.81317.41414.91712.42910.996A=∑At024.49633.24937.38139.70741.88742.69244.496Speedup1.0001.3841.5091.6511.7151.7451.7901.789Threads910111213141516At=∑Ak9.8548.8617.9627.2966.6516.0825.8035.486A=∑At45.82945.24944.60845.31745.63245.67446.64947.515Speedup1.7801.8201.8821.8861.8901.9081.8881.868


## 4.3. Application of Constraints

In this work, only the Dirichlet constraint is considered. The parallel algorithm for Dirichlet constraints has two parts as shown in Figure5. The first part is to remove all the Dirichlet rows from A while setting its diagonal entry to 1. In this experiment, a two-level nested std:map is used to store a sparse matrix. The std:map is not a thread-safe data structure. The removal of rows means the removal of internal std:map object from outer std:map object where each edge represents a row. The second part is more time-consuming as for every non-Dirichlet row, Dirichlet columns entries are searched and removed. Figure 6 shows a 7 × 7 stiffness matrix before and after applying Dirichlet constraints where nodes 4 and 6 lay on the Dirichlet boundary. Removal of Dirichlet rows was performed sequentially because of a non-thread-safe data structure. Non-Diriechlet rows were distributed among the threads where each thread will remove Dirichlet columns from each row. Table 10 shows timing and speedup for the Dirichlet algorithm. The second part of the algorithm is embarrassingly parallel and has achieved an almost linear speedup of 14.65%for 16 threads.Figure 5
Parallel Dirichlet implementation.Figure 6
Dirichlet constraints. (a) Global stiffness matrix. (b) Global stiffness matrix after applying Dirichlet constraints.
(a)(b)Table 10
Dirichlet application using multiple threads using a problem having 2099201 unknown.
Threads12345678Time1671825551423343281241218Speedup1.002.033.033.954.875.966.927.68Threads910111213141516Time194173158149138128119114Speedup8.619.6810.5811.2512.1613.1114.0614.65


## 4.4. Linear Algebra Solvers

There were two most commonly used iterative solvers Jacobi and Conjugate Gradient considered in this work. These are implemented in parallel using OpenMP. Their algorithms and implementation details are as follows. All the critical sections are colored blue while synchronization points are colored red.(1)
Jacobi solver iteration can be segregated into the following parts.(i)
Calculation of residualr=b−Ax(ii)
New approximation ofx+=D−1r(iii)
Calculation of error thresholde=r/x(2)
Conjugate Gradient SolverAlgorithm 1: Jacobi solver.
Input:A, b⟶,x, N, x_norm, r_normOutput:x(1)
fork (A, b⟶,x, N, x_norm, r_norm)(2)
Th_id, Tt_id (Thread id and Total Thread by OS)(3)
St_RowId ← Ed_RowId ← 0 {Start and end ids}(4)
IfTh_id ＜ N mod Tt_ththen(5)
St_RowID ← Th_id∗N/Tt/th+Th_id(6)
Ed_RowID ← St_RowID+N/Tt_th+1(7)
else(8)
St_RowID ← Th_id ∗ N/Tt_th+Tt_th(9)
Ed_RowID ← St_RowID+N/Tt_th(10)
end_if(11)
r {Local Residual vector for allocated nodes}(12)
error ← 1(13)
Whileerror > Thresholddo(14)
r_tmp ← x_tmp ← 0 {Local Norms for thread}(15)
foreachi ← St_RowID, Ed,RowIDdo(16)
repeat(17)
r[i-St_RowID]←b[i]-A[i,j]∗x[j](18)
untilFor each entry j for row i of A(19)
r_tmp ← r_tmp + r[i-St_RowID]2(20)
end for(21)
x_norm ← r_norm ← 0 [Global Norms](22)
barrier(23)
foreachi ← St_RowID, Ed_RowIDdo(24)
x[i] ← x[i] + r[i-St_RowID]/A[i,i](25)
x_tmp ← x_tmp+x[i]2(26)
end for(27)
lock(28)
r_norm ←r_norm+r_tmp(29)
x_norm←x_norm+x_tmp(30)
unlock(31)
barrier(32)
error ←r_norm/x_norm(33)
end while(34)
joinFigure 7
Parallel Jacobi solver stages. (a) Parallel residual calculation. (b) Error calculation.The second implementation for Jacobi is to create a single parallel region that will cover all three parts of Jacobi iteration as shown in Algorithm1. There is only one parallel region in this application starting at line 1 and collapsing at line 35. All three stages are implemented in single parallel regions, and different stages have to be segregated using explicit barriers at line 22. The first explicit barrier is necessary as in the first region x is used to calculate r and in the second region x is updated using r. The second explicit barrier at line 31 is required to the calculation of error which requires r/x. Before the second barrier, a critical section, from lines 27 to 30, is required to calculate r and x. Table 11 shows timing data for Algorithm 1 calculated using a mesh of 131585 nodes and 263168 elements. Jacobi has a slow convergence rate that is taking a large number of iterations to converge and it is the most time-consuming part of FE solver. For this problem, it took about 165009 iterations. The parallel task has enabled the reduction of the Jacobi solver's sequential time of 3780 sec to parallel time of 605 sec using 16 threads and achieving 6.246 times speedup. During execution, the stopping criteria adopted was either the error is less than 10–8 or total iterations are greater than 106.Table 11
Jacobi solver parallel timing and speedup for system of equations having 131585 unknowns.
Threads12345678Jacobi3780258320231458145511781129982Speedup1.001.461.872.592.603.213.353.85Threads910111213141516Jacobi899814754733681642619605Speedup4.204.645.015.165.555.896.116.25Algorithm2 shows a parallel algorithm of Conjugate Gradient. While designing this algorithm, the main aim was to keep the global synchronization minimum. In this algorithm, global synchronization is represented by barrier function at lines 17, 25, and 34. All the codes displayed in the black color are each thread individual code that has to perform without communicating with other threads. All the codes in blue are a critical section that has to be executed by all threads but one at a time. In this section, global variables rtr1, rtr2, and ptAp are computed. In each iteration, there are two synchronization points required to update shared variables such as ptAp and rtr2, respectively.Algorithm 2: Parallel Conjugate Gradient algorithm.
Input: (A, b⟶,x, N, rtt1, rtr2, ptAp)Output:x(1)
fork (A, b⟶,x, N, rtt1, rtr2, ptAp) {shared variables}(2)
Th_id, Tt_id (Thread id and Total Thread by OS)(3)
St_RowId ← Ed_RowId ← 0 {Start and end ids}(4)
IfTh_id ＜ N % Tt_ththen(5)
St_RowID ← Th_id∗N/Tt+Th_id(6)
Ed_RowID ← St_RowID+N/Tt_th+1(7)
else(8)
St_RowID ← Th_id ∗ N/Tt_th+Tt_th(9)
Ed_RowID ← St_RowID+N/Tt_th(10)
end_if(11)
rtr1 ← rtr2 ← t_rtr ← 0(12)
foreachi ← St_RowID.Ed_RowIDdo(13)
p[i] ← r[i] ← b[i]-Ai:,x(14)
t_rtr ← t_rtr + r[i] ∗ r[i](15)
end for(16)
lockrtr2 ← rtr2 + t_rtr unlock(17)
barrier(18)
whilertr2>Threshold do(19)
rtr1 ← rtr2(20)
rtr2 ← t_ptAp ← 0(21)
foreachi ← St_RowID, Ed_RowIDdo(22)
t_p'Ap ← t_ptAp+p[i].Ai:,p(23)
end for(24)
lockptAp ← ptAp+t_ptAp unlock(25)
barrier(26)
t_rtr ← 0(27)
foreachi ← St_RowID, Ed_RowIDdo(28)
x[i] ← x[i]+(rtr1/ptAp).p[i](29)
r[i] ← r[i]+(rtr1/ptAp).Ai:,p(30)
t_rtr ← t_rtr+r[i] ∗ r [i](31)
end for(32)
ptAp ← 0(33)
lockrtr2 ← rtr2+t_rtr unlock(34)
barrier(35)
foreachi ← St_RowID, Ed_RowIDdo(36)
p[i] ← r[i]+(rtr2/rtr1).p[i](37)
end for(38)
end while(39)
joinTables12 and 13 show the timing and speedup data of Algorithm 2 for two mashes. The system of 131585 unknown took 1019 iterations while the system of 2099281 unknown took 3987 iterations. Algorithm 2 has three critical sections (all colored in blue) where the first critical section at line 16 has to be executed once as it is before while loop while the other two sections at lines 24 and 33 are in a while loop and have to be executed in every iteration. These sections are to update global variables (rtr and ptAp) with thread-local calculations. These variables cannot be updated in the same critical section as ptAp to use for the calculation of rtr. Therefore, every time a thread after going through a critical section has to wait for other threads to execute it before further processing, this is imposed by global communication. For a fixed size problem, as the number of threads in the parallel region increases, the time required for thread-independent tasks (colored in black grey) will reduce while the critical section and communication timing will increase.Table 12
Conjugate Gradient solver parallel timing and speedup for system of equations having 131585 unknowns.
Threads12345678CG23.4217.2612.6410.108.437.296.906.10Speedup1.001.361.852.322.783.213.393.85Threads910111213141516CG5.465.124.694.394.233.853.653.62Speedup4.294.575.005.345.546.096.426.47Table 13
Conjugate Gradient solver parallel timing and speedup for system of equations having 2099201 unknowns.
Threads12345678CG15661214997931721687670634Speedup1.001.291.571.682.172.282.342.47Threads910111213141516CG520519.95479449399366.93366312Speedup3.013.013.273.493.934.274.285.01


## 5. Parallel FE Solver Performance

The overall application performance obtained during experiments is given in Tables14 and 15. These tables show the time consumed to solve the finite element problem with 131585 and 2099201 unknowns.Table 14
Finite element solver timing, speedup, and efficiency for system of 131585 unknowns.
Threads12345678Time26.2319.0113.9111.059.338.157.756.94Speedup1.001.381.892.372.813.223.393.78Efficiency1.000.690.630.590.560.540.480.47Threads910111213141516Time6.295.965.525.205.074.694.514.50Speedup4.174.404.755.045.185.595.825.83Efficiency0.460.440.430.420.400.400.390.36Table 15
Finite element solver timing, speedup, and efficiency for system of 2099201 unknowns.
Threads12345678Time167112661031956743706686649Speedup1.001.321.621.752.252.372.4352.58Efficiency1.000.660.540.440.450.400.3480.32Threads910111213141516Time534532490460408377374321Speedup3.133.143.413.644.094.444.4685.21Efficiency0.350.310.310.300.320.320.2980.33Figure8 shows speedup for three problems (smaller to larger from left to right) using threads 1 to 16 using Conjugate solver. The data structure used to collect nodes and elements is not thread-safe. It can be seen that, for smaller mash, speedup has started to decline for more threads while as the number of unknowns increases, speedup keeps on increasing for higher threads as well. Construction of equation and Dirichlet speedup has almost linearly increased with an increase in threads. In both of these cases, the choice of data container for A has a significant impact as discussed previously. The last stage of the FE solver has taken the maximum time and achieved the least speedup among all. Iterative solvers (Jacobi and Conjugate Gradient) have a major impact on solver complete time, so the increase in solver speedup is similar to iterative speedup.Figure 8
FE solver of all stages’ speedup achieved for 16 threads. There are three different size problems displayed.
(a)(b)(c)


## 6. Conclusion

This work has focused on the implementation of Finite Element Application (FEA) using object-oriented programming for a shared memory architecture using OpenMP without considering conventional techniques of application speedup and specialized data structures like CSR, etc. FEA has been segregated into four stages which are read mesh, construction of equations, application of boundary conditions, and solution computation for a linear system of equations. The models available for solving FEM-associated issues work very similarly to each other. The key difference among these models lies in the fact that every model uses a different approach to handle the complicated geometries (and boundaries) with relative case All timings are calculated on a computer having Xeon E5-2650 (octa-core and support hyperthreading) processor and 32 GB RAM. In the initial stage, the mesh is read, in the most efficient block of nodes, and elements (block partitioning) are read in each thread using a local file descriptor and are stored in a shared data container. In the second stage, a system of equations is constructed from nodes and elements. Standard template library std:map class is used for sparse matrix container which is not a thread-safe data structure. Sparse matrix container is implemented using a two-level std:map where an individual std:map object is used for each row and an individual lock is associated with each row. If all threads update different rows, then no thread will have to wait as each row has a separate object and lock while on the other extreme, if all threads write to the same row, only one will have access to the row object while others have to wait. This stage has a maximum speedup of 12.67 times. In the next stage, the Dirichlet condition is implemented in parallel and achieved a maximum speedup of 14.65 times. In the last stage, two linear algebra solvers (Jacobi and Conjugate Gradient) are implemented to compute the solution. Jacobi has two synchronization points per iteration while Conjugate Gradient has three. Conjugate Gradient has a higher level of complexity and a fast divergence rate as compared to Jacobi. Parallel Jacobi has achieved 6 times speedup to its sequential Jacobi implementation. Parallel Conjugate Gradient has achieved 6 times speedup to its sequential Conjugate Gradient. The last stage is the most time-consuming part of FEA and overall speedup directly depends on this stage. During experiments, it has been observed that about 85% to 90% of the total time was consumed by solvers. To the best of our knowledge, in the domain of shared memory architecture, solving the issues related to FEM while designing multithreaded architecture is the pioneering work and this work still has potential for improvement.



---

*Source: 2375102-2022-03-19.xml*