Academia.eduAcademia.edu

Figure 13 – uploaded by Mohamed H. Aissa

See full PDFdownloadDownload figure
Table 6.3: Execution Time of the eDisco on GPU and different number of CPU cores combination

Table 6 3: Execution Time of the eDisco on GPU and different number of CPU cores combination

Related Figures (62)

Figure 2.1: Generic Diagram of the Single Instruction Single Data System (SISD) This category includes single-processor computers based on the common Von Neumann model illustrated in Figure 2.1. This is a serial processing hardware architecture.
Figure 2.1: Generic Diagram of the Single Instruction Single Data System (SISD) This category includes single-processor computers based on the common Von Neumann model illustrated in Figure 2.1. This is a serial processing hardware architecture.
Figure 2.2: Generic Diagram of the Single Instruction Multiple Data System (SIMD), based on [1]
Figure 2.2: Generic Diagram of the Single Instruction Multiple Data System (SIMD), based on [1]
Figure 2.3: Generic Diagram of the Multiple Instructions Multiple Data (MIMD),based on [1]
Figure 2.3: Generic Diagram of the Multiple Instructions Multiple Data (MIMD),based on [1]
The first class of shared-memory architecture (UMA) is also called symmetric multi- processor (SMP)/1][15]. The access of all memory locations occurs at the same speed, which causes the architecture to be symetric. This is indeed the simplest shared-memory architecture since all processors share the same connection (bus system) to the shared memory as shown in Figure 2.4. The bus system can not serve a large number of proces- sors, unless the memory access speed declines dramatically. The memory bandwidth is a limiting factor for UMA architecture. Figure 2.4: Generic Diagram of the Uniform Memory Access (UMA) Architecture
The first class of shared-memory architecture (UMA) is also called symmetric multi- processor (SMP)/1][15]. The access of all memory locations occurs at the same speed, which causes the architecture to be symetric. This is indeed the simplest shared-memory architecture since all processors share the same connection (bus system) to the shared memory as shown in Figure 2.4. The bus system can not serve a large number of proces- sors, unless the memory access speed declines dramatically. The memory bandwidth is a limiting factor for UMA architecture. Figure 2.4: Generic Diagram of the Uniform Memory Access (UMA) Architecture
Figure 2.5: Generic Diagram of the Non-Uniform Memory Access (NUMA) Architecture The non-uniform access is reflected in the access speed of the shared memory. Physically near shared-memory parts are accessed faster then other parts. With a higher number of processors the complexity of the NUMA architecture increases leading to important memory access delays known as latency. One strategy to minimize the memory traffic is to use a cache associated with every processor. The cache entries are kept coherent! through a protocol. This modified architecture is called cache coherent NUMA or ccNUMA. This modification leads however to a more complicated bus structure which raises the costs [15].
Figure 2.5: Generic Diagram of the Non-Uniform Memory Access (NUMA) Architecture The non-uniform access is reflected in the access speed of the shared memory. Physically near shared-memory parts are accessed faster then other parts. With a higher number of processors the complexity of the NUMA architecture increases leading to important memory access delays known as latency. One strategy to minimize the memory traffic is to use a cache associated with every processor. The cache entries are kept coherent! through a protocol. This modified architecture is called cache coherent NUMA or ccNUMA. This modification leads however to a more complicated bus structure which raises the costs [15].
Figure 2.6: Generic Diagram of the Distributed-Memory Architecture 2.1.3 Hybrid Architecture
Figure 2.6: Generic Diagram of the Distributed-Memory Architecture 2.1.3 Hybrid Architecture
Figure 2.7: Generic Diagram of the Distributed-Shared Memory Architecture
Figure 2.7: Generic Diagram of the Distributed-Shared Memory Architecture
Figure 2.8: Generic Diagram of the GPU Cluster Architecture NVIDIA, the other major player in the GPU solution for HPC, has developed three architectures: Tesla, Fermi and Kepler (see table 2.1). The latter architecture was originally planned to be released late in 2011. However the release has been postponed until 2012 [15]. Like AMD, NVIDIA provides an SDK with CUDA as programming language (see 3.2).
Figure 2.8: Generic Diagram of the GPU Cluster Architecture NVIDIA, the other major player in the GPU solution for HPC, has developed three architectures: Tesla, Fermi and Kepler (see table 2.1). The latter architecture was originally planned to be released late in 2011. However the release has been postponed until 2012 [15]. Like AMD, NVIDIA provides an SDK with CUDA as programming language (see 3.2).
Figure 2.9: Summary Diagram of mentioned Architectures All the mentioned architectures are depicted in Figure 2.9. Flynn’s taxonomy was used to map the different architectures based on the data and instruction streams. The MIMD category was further split according to the memory organization into shared-memory, distributed memory and distributed shared-memory. The computational accelerators has been described with focus on the graphic processing unit (GPU). Colored ones represents used architecture in this work: uniform access memory system and graphic processing unit (see the implementation chapter 5).
Figure 2.9: Summary Diagram of mentioned Architectures All the mentioned architectures are depicted in Figure 2.9. Flynn’s taxonomy was used to map the different architectures based on the data and instruction streams. The MIMD category was further split according to the memory organization into shared-memory, distributed memory and distributed shared-memory. The computational accelerators has been described with focus on the graphic processing unit (GPU). Colored ones represents used architecture in this work: uniform access memory system and graphic processing unit (see the implementation chapter 5).
Figure 2.10: OpenMP fork-join execution model The OpenMP uses the fork-join model of parallel execution [20]. As illustrated in Figure 2.10 OpenMP program begins as a single thread of execution, called the initial thread. When any thread encounters a parallel region marked by the compiler directives, the thread forks additional threads creating a team of threads and becomes the master of the new team. The master thread and the new created threads execute the desired program part in parallel. There is an implicit synchronization barrier at the end of the parallel region. The threads wait until the full team reaches the barrier, and then they join back together. At that point the original or master thread continues until the next parallel region or the end of the program. The clause nowait deactivates the implicit synchronization through the barrier. This barrier forces all threads to wait for the slowest one among them.
Figure 2.10: OpenMP fork-join execution model The OpenMP uses the fork-join model of parallel execution [20]. As illustrated in Figure 2.10 OpenMP program begins as a single thread of execution, called the initial thread. When any thread encounters a parallel region marked by the compiler directives, the thread forks additional threads creating a team of threads and becomes the master of the new team. The master thread and the new created threads execute the desired program part in parallel. There is an implicit synchronization barrier at the end of the parallel region. The threads wait until the full team reaches the barrier, and then they join back together. At that point the original or master thread continues until the next parallel region or the end of the program. The clause nowait deactivates the implicit synchronization through the barrier. This barrier forces all threads to wait for the slowest one among them.
Figure 2.11: Generic Diagram of the Cluster of Nodes Architecture with MPI and OpenMP Models
Figure 2.11: Generic Diagram of the Cluster of Nodes Architecture with MPI and OpenMP Models
Figure 2.12: Diagram of the parallel programming design stages (based on [:
Figure 2.12: Diagram of the parallel programming design stages (based on [:
Figure 2.13: Diagram of different speedups and processors numbers After all cells have calculated their temperatures, boundary area (blue and mauve right on Figure 2.13) has to be updated. This update is performed on a large number of cells. The Figure 2.13 illustrates a periodic boundary condition where the state at the end of the computation domain has to be copied at the end of the domain.
Figure 2.13: Diagram of different speedups and processors numbers After all cells have calculated their temperatures, boundary area (blue and mauve right on Figure 2.13) has to be updated. This update is performed on a large number of cells. The Figure 2.13 illustrates a periodic boundary condition where the state at the end of the computation domain has to be copied at the end of the domain.
Figure 2.15: Diagram of different speedups and processors numbers A program comprise a serial portion and a portion that could be parallelized. A parallel portion runs two times faster on two-processor systems as in one-processor system. The more processors are used the faster the parallel part will run. The serial part execution time remains untouched since it executes only on one processor. The Figure 2.15 illustrates the execution time of a program running on different numbers of processors. One fifth of the program is executed serially. The execution time of the serial part is constant and remains one fifth of the total execution time. The parallel execution time is inversely proportional to the number of processors. Four processors for instance are four time faster in running the parallel part but only two and half time faster than one processor running the entire program containing a serial portion of one fifth. The amount of performance benefit an application will realize by running on parallel systems depends entirely on the extent to which it can be parallelized [22]. G. Amdahl [23] was the first to model the performance expectation as follow :
Figure 2.15: Diagram of different speedups and processors numbers A program comprise a serial portion and a portion that could be parallelized. A parallel portion runs two times faster on two-processor systems as in one-processor system. The more processors are used the faster the parallel part will run. The serial part execution time remains untouched since it executes only on one processor. The Figure 2.15 illustrates the execution time of a program running on different numbers of processors. One fifth of the program is executed serially. The execution time of the serial part is constant and remains one fifth of the total execution time. The parallel execution time is inversely proportional to the number of processors. Four processors for instance are four time faster in running the parallel part but only two and half time faster than one processor running the entire program containing a serial portion of one fifth. The amount of performance benefit an application will realize by running on parallel systems depends entirely on the extent to which it can be parallelized [22]. G. Amdahl [23] was the first to model the performance expectation as follow :
Figure 2.16: Expected speedup as a function of the parallel portion of a program running on 100 processors
Figure 2.16: Expected speedup as a function of the parallel portion of a program running on 100 processors
Figure 2.17: Expected speed up as a function of the number of processors for different parallel potions
Figure 2.17: Expected speed up as a function of the number of processors for different parallel potions
Figure 3.1: Comparison of the GPU and CPU architecture [3] The architecture of graphics computing units differs from the CPU architecture as the GPU was originally designed for graphics generation. Unlike the CPU most of the capacities are used for the computing. The CPU performs a multitude of different tasks, therefore most of the capabilities are invested on the memory caching and the control flow as shown in Figure 3.1. The GPU however contains more arithmetic logical units (ALU) dedicated for floating-point calculation. Graphical processing units are therefore taking more advantage from the Moore’s Law predicting the growth of the number of transistors on a chip.
Figure 3.1: Comparison of the GPU and CPU architecture [3] The architecture of graphics computing units differs from the CPU architecture as the GPU was originally designed for graphics generation. Unlike the CPU most of the capacities are used for the computing. The CPU performs a multitude of different tasks, therefore most of the capabilities are invested on the memory caching and the control flow as shown in Figure 3.1. The GPU however contains more arithmetic logical units (ALU) dedicated for floating-point calculation. Graphical processing units are therefore taking more advantage from the Moore’s Law predicting the growth of the number of transistors on a chip.
Block consists of 1 to 512 concurrent threads for GPUs with Tesla architecture [3]. Each thread has a unique thread ID (TID) starting with 0. Figure 3.2: illustration of the GPU device architecture [4]
Block consists of 1 to 512 concurrent threads for GPUs with Tesla architecture [3]. Each thread has a unique thread ID (TID) starting with 0. Figure 3.2: illustration of the GPU device architecture [4]
from all others. The order of running the blocks is therefore not relevant. If two blocks or more blocks executes in parallel the result remains the same. The programmer declares the number of threads of a block. The resource demand in form of registers and shared memory depends on the code being executed within a thread. When a streaming multiprocessor has enough resources for the block, the block of threads is created and every thread receives a thread ID (TID). Otherwise the block fails to start and the user is informed through a run-time error message (see Subsection 3.2.2). A streaming multiprocessor can start therefore at least one block of threads and up to eight concurrently [6]. The Tesla C1060 can thus execute from eight to sixty four blocks concurrently.
from all others. The order of running the blocks is therefore not relevant. If two blocks or more blocks executes in parallel the result remains the same. The programmer declares the number of threads of a block. The resource demand in form of registers and shared memory depends on the code being executed within a thread. When a streaming multiprocessor has enough resources for the block, the block of threads is created and every thread receives a thread ID (TID). Otherwise the block fails to start and the user is informed through a run-time error message (see Subsection 3.2.2). A streaming multiprocessor can start therefore at least one block of threads and up to eight concurrently [6]. The Tesla C1060 can thus execute from eight to sixty four blocks concurrently.
Figure 3.4: Illustration of SIMT Warps Scheduling of two Blocks of Threads A SIMT streaming multiprocessor realizes full efficiency and performance when all 32 threads of a warp follow the same execution path [6]. If threads of a warp diverge via a data-dependent conditional branch, the warp serially executes each branch path taken as illustrated in Figure 3.5, disabling threads that are not on that path. When all paths complete, the threads reconverge to the original execution path. Thread divergence affects only threads within the same warp. Different warp can execute different code paths.
Figure 3.4: Illustration of SIMT Warps Scheduling of two Blocks of Threads A SIMT streaming multiprocessor realizes full efficiency and performance when all 32 threads of a warp follow the same execution path [6]. If threads of a warp diverge via a data-dependent conditional branch, the warp serially executes each branch path taken as illustrated in Figure 3.5, disabling threads that are not on that path. When all paths complete, the threads reconverge to the original execution path. Thread divergence affects only threads within the same warp. Different warp can execute different code paths.
Figure 3.5: Illustration of the serial execution of a divergent warp of threads For program correctness, programmers can essentially ignore SIMT execution attributes such as warps. They can achieve however substantial performance improvements by writing code that seldom requires threads in a warp to diverge.
Figure 3.5: Illustration of the serial execution of a divergent warp of threads For program correctness, programmers can essentially ignore SIMT execution attributes such as warps. They can achieve however substantial performance improvements by writing code that seldom requires threads in a warp to diverge.
on threads. ] Every thread runs separately a copy of the kernel. CUDA provides built-in variables? lis ted in table 3.1, which combined generate a thread identifier (tid) unique to every launched thread (line 3 Listing 3.2). The conditional construct surrounding the add operatio perform the n (line5-7) ensure through the condition (tid<N) that only needed threads add operation?. Threads are created with the kernel launch from the host. Listing 3.3 shows t call of the add function. It differs semantically from C/C++ function calls use of new execution configuration syntax (<<<,>>>). Between the inequa the user gives the number of blocks and the number of threads that should be These given parameters are saved then automatically in the built-in variables he kernel with the ity signs aunched gridDim and blockDim as listed in table 3.1. The kernel launch is asynchronous, which means the host can proceed other tasks on CPU directly after the kernel call before t computation is finished. he kernel
on threads. ] Every thread runs separately a copy of the kernel. CUDA provides built-in variables? lis ted in table 3.1, which combined generate a thread identifier (tid) unique to every launched thread (line 3 Listing 3.2). The conditional construct surrounding the add operatio perform the n (line5-7) ensure through the condition (tid<N) that only needed threads add operation?. Threads are created with the kernel launch from the host. Listing 3.3 shows t call of the add function. It differs semantically from C/C++ function calls use of new execution configuration syntax (<<<,>>>). Between the inequa the user gives the number of blocks and the number of threads that should be These given parameters are saved then automatically in the built-in variables he kernel with the ity signs aunched gridDim and blockDim as listed in table 3.1. The kernel launch is asynchronous, which means the host can proceed other tasks on CPU directly after the kernel call before t computation is finished. he kernel
kernel can call a device function, which run on the GPU device and can be called only within the device from a kernel or another device function. The declaration of this type of functions starts with a qualifier “___ device____” (See table 3.2). Kernels are executed in the device but callable only from the host. A kernel cannot return a value. A device function cannot call a kernel. A host function could be the main function or any other ANSI C conventional function running on the CPU host.
kernel can call a device function, which run on the GPU device and can be called only within the device from a kernel or another device function. The declaration of this type of functions starts with a qualifier “___ device____” (See table 3.2). Kernels are executed in the device but callable only from the host. A kernel cannot return a value. A device function cannot call a kernel. A host function could be the main function or any other ANSI C conventional function running on the CPU host.
Figure 3.6: Overview of the CUDA device memory model|5] 3 Use of Graphics processing Units for Parallel Computation The API function cudaMalloc() allocates systematically the needed memory in the off-chir global memory. This type of memory is indeed 4GB large (see table A.1) and visible to all threads of all blocks but offers a slow access (latency). A read/write operation lasts between 400 and 800 clock cycles[3]. A variable lying in the global memory can be read and modified from any kernel or from the host by some API functions like cudaMemcpy The content of those variables exists during the whole application and remains untouched after one kernel is finished.
Figure 3.6: Overview of the CUDA device memory model|5] 3 Use of Graphics processing Units for Parallel Computation The API function cudaMalloc() allocates systematically the needed memory in the off-chir global memory. This type of memory is indeed 4GB large (see table A.1) and visible to all threads of all blocks but offers a slow access (latency). A read/write operation lasts between 400 and 800 clock cycles[3]. A variable lying in the global memory can be read and modified from any kernel or from the host by some API functions like cudaMemcpy The content of those variables exists during the whole application and remains untouched after one kernel is finished.
To encounter the latency problem CUDA provides other type of memories that are not so large as the global memory and not overall visible, therefore faster accessible (see table 3.3).
To encounter the latency problem CUDA provides other type of memories that are not so large as the global memory and not overall visible, therefore faster accessible (see table 3.3).
Figure 4.1: One-dimensional Illustration of the inter-elements Discontinuity of the DG Global Solution The local DG approach of element-based basis functions leads to a discontinuous approximation solution q,. In each element interface there are two different solutions, one emerging from the element itself (q,,) and the other from the neighbor element (qp,_,, ). An illustration of the one-dimensional discontinuity of the local solutions q;,, is depicted in Figure 4.1.
Figure 4.1: One-dimensional Illustration of the inter-elements Discontinuity of the DG Global Solution The local DG approach of element-based basis functions leads to a discontinuous approximation solution q,. In each element interface there are two different solutions, one emerging from the element itself (q,,) and the other from the neighbor element (qp,_,, ). An illustration of the one-dimensional discontinuity of the local solutions q;,, is depicted in Figure 4.1.
Table 4.1: RK Time Discretization Parameters where aj; and 3; are in Table 4.1 defined [35] Starting from the above described method a three-stage three-order exact Runge-Kutta method is used to integrate (4.13),cf. [29]:
Table 4.1: RK Time Discretization Parameters where aj; and 3; are in Table 4.1 defined [35] Starting from the above described method a three-stage three-order exact Runge-Kutta method is used to integrate (4.13),cf. [29]:
Figure 4.2: Illustration of a numerical domain subdivided in interior, Boundary and Ghost Elements
Figure 4.2: Illustration of a numerical domain subdivided in interior, Boundary and Ghost Elements
The solver routine is a function that uses the initial solution and the grid to solve the governing equations numerically and to produce a solution for a given time point t = tfinal. The solver, as illustrated in Algorithm 4.1, executes a series of iterations on the governing equations to proceed in time until the final time is reached. In every iteration a Runge-Kutta multiple step scheme (lines 3-6) is executed to solve the equations for one time step. The RK scheme requires a residual computation to update the physical solution (line 4). The eDisco program computes the residual for every time step using a set of functions solving the equation (4.9), as shown in Algorithm 4.2. One function computes the first term based on a surface integral S$? (line 1) and the second function adds a second term based on a volume integral H}. (line 2). The result is then multiplied with the inverse mass matrix M, i ,, (line 3). The final solution reached after a fixed amount of time stepping is stored in an output file.
The solver routine is a function that uses the initial solution and the grid to solve the governing equations numerically and to produce a solution for a given time point t = tfinal. The solver, as illustrated in Algorithm 4.1, executes a series of iterations on the governing equations to proceed in time until the final time is reached. In every iteration a Runge-Kutta multiple step scheme (lines 3-6) is executed to solve the equations for one time step. The RK scheme requires a residual computation to update the physical solution (line 4). The eDisco program computes the residual for every time step using a set of functions solving the equation (4.9), as shown in Algorithm 4.2. One function computes the first term based on a surface integral S$? (line 1) and the second function adds a second term based on a volume integral H}. (line 2). The result is then multiplied with the inverse mass matrix M, i ,, (line 3). The final solution reached after a fixed amount of time stepping is stored in an output file.
The surface integral and the volume integral are essential components of the DG solver and therefore described in the following Algorithms (4.3,4.4). The volume integral is a loop over all elements of the computational domain (lines 1-10). For every element the contributions of all integration points are calculated (lines 3-9) and added to the integral value of the element in every degree of freedom (line 9).
The surface integral and the volume integral are essential components of the DG solver and therefore described in the following Algorithms (4.3,4.4). The volume integral is a loop over all elements of the computational domain (lines 1-10). For every element the contributions of all integration points are calculated (lines 3-9) and added to the integral value of the element in every degree of freedom (line 9).
The surface integral on the other hand runs over all the edges of the computational domain (line 1-9). For every edge the contributions of the integration points are calculated and added to the related degree of freedom of the element in the right (line 9) and the left (line 8) of the edge.
The surface integral on the other hand runs over all the edges of the computational domain (line 1-9). For every edge the contributions of the integration points are calculated and added to the related degree of freedom of the element in the right (line 9) and the left (line 8) of the edge.
Algorithm 4.4 Algorithm of the Volume Integral of The DG approach
Algorithm 4.4 Algorithm of the Volume Integral of The DG approach
Figure 5.1: Illustration of the four Steps of the Design Cycle for GPU Applications
Figure 5.1: Illustration of the four Steps of the Design Cycle for GPU Applications
Figure 5.3: Illustration of the parallelization strategy
Figure 5.3: Illustration of the parallelization strategy
Figure 5.4: Illustration of the parallelization of two function separated by a function running o: CPU
Figure 5.4: Illustration of the parallelization of two function separated by a function running o: CPU
Listing 5.1: Implementing of a CPU for-loop (left) on CUDA (right):
Listing 5.1: Implementing of a CPU for-loop (left) on CUDA (right):
As important as the thread management the kernel itself has to be considered. A kerne could be fully parallel as illustrated in Figure 5.6. The kernel starts an amount! o threads, which run independently performing parallel computation and at the end save the results in separate memory positions. No memory conflict occurs. Some functions are however not fully parallel. The volume integral computes fluxes on integration points and save the results in every degree of freedom. This leads to thread conflicts since different threads compute different fluxes on different integration points and save their result in Nested loops are also a common parallelization pattern. The easiest way is to start two-dimensional blocks of threads as seen in Algorithm 5.1. The two indices blockIdx.x and blockIdx.y replace the two loops indexes i and j. The code could remain the same if no data dependency is identified. Synchronization may be required depending on the task. In comparison with the CPU algorithm 4.3 of the volume integral, the GPU algorithm replaces the two first loops covering the elements and the integration points. A third FOR-loop covering the degrees of freedom has been kept because of the high data dependency among the DOFs.
As important as the thread management the kernel itself has to be considered. A kerne could be fully parallel as illustrated in Figure 5.6. The kernel starts an amount! o threads, which run independently performing parallel computation and at the end save the results in separate memory positions. No memory conflict occurs. Some functions are however not fully parallel. The volume integral computes fluxes on integration points and save the results in every degree of freedom. This leads to thread conflicts since different threads compute different fluxes on different integration points and save their result in Nested loops are also a common parallelization pattern. The easiest way is to start two-dimensional blocks of threads as seen in Algorithm 5.1. The two indices blockIdx.x and blockIdx.y replace the two loops indexes i and j. The code could remain the same if no data dependency is identified. Synchronization may be required depending on the task. In comparison with the CPU algorithm 4.3 of the volume integral, the GPU algorithm replaces the two first loops covering the elements and the integration points. A third FOR-loop covering the degrees of freedom has been kept because of the high data dependency among the DOFs.
Figure 5.5: Degree of freedom and integration Points for Surface and Volume Integral on a Triangle Element with first polynomial order (p=1) the same degree of freedom. A first polynomial-order triangle element is illustrated in Figure 5.5. For the volume integral three integration points situated at the center of every edge are used to compute the integral on every degree of freedom. The contributions of different integration points are computed in parallel. The contri- butions are added at the end of the integral and saved in the degrees of freedom. As shown in Figure 5.6 the kernel starts a set of threads, one per integration point. Every three threads are computing the integral on the same degree of freedom. The addition has to be performed serially otherwise a race condition emerges. Refer to Section 5.5 for the impact of this error on the final result. A full parallel kernel can be started after a partial parallel kernel as every kernel gets its new execution context? and parameter while starting.
Figure 5.5: Degree of freedom and integration Points for Surface and Volume Integral on a Triangle Element with first polynomial order (p=1) the same degree of freedom. A first polynomial-order triangle element is illustrated in Figure 5.5. For the volume integral three integration points situated at the center of every edge are used to compute the integral on every degree of freedom. The contributions of different integration points are computed in parallel. The contri- butions are added at the end of the integral and saved in the degrees of freedom. As shown in Figure 5.6 the kernel starts a set of threads, one per integration point. Every three threads are computing the integral on the same degree of freedom. The addition has to be performed serially otherwise a race condition emerges. Refer to Section 5.5 for the impact of this error on the final result. A full parallel kernel can be started after a partial parallel kernel as every kernel gets its new execution context? and parameter while starting.
Figure 5.6: Illustration of the execution of full parallel (left) and partial parallel kernels (right
Figure 5.6: Illustration of the execution of full parallel (left) and partial parallel kernels (right
Every element is an array of degrees of freedom and every degree of freedom is a structure with the conservative flow variables (p, pu, pv,e). For the graphics computational units the data has been rearranged in one-dimensional arrays. The access of the data is then performed using array index manipulation:
Every element is an array of degrees of freedom and every degree of freedom is a structure with the conservative flow variables (p, pu, pv,e). For the graphics computational units the data has been rearranged in one-dimensional arrays. The access of the data is then performed using array index manipulation:
5.4 Implemented Functions
5.4 Implemented Functions
Figure 5.9: Flux calculation on Element Edges The race condition problem can be also hidden in the algorithm. The surface integral in this work, for instance, was implemented in the CPU in a way causing a hard-to-find race condition. The integral is calculated for every edge and added to the two neighbor elements on the edge as seen in Figure 5.9. Because the fluxes are computed in parallel, a case occurs that two or three different edges are writing in the same elements.
Figure 5.9: Flux calculation on Element Edges The race condition problem can be also hidden in the algorithm. The surface integral in this work, for instance, was implemented in the CPU in a way causing a hard-to-find race condition. The integral is calculated for every edge and added to the two neighbor elements on the edge as seen in Figure 5.9. Because the fluxes are computed in parallel, a case occurs that two or three different edges are writing in the same elements.
Figure 5.10: Effect of the Race Condition on an acoustic pulse results: A race condition free output is shown on the left figure. The right figure is the output of an implementation with a race condition causing a noise in the results.
Figure 5.10: Effect of the Race Condition on an acoustic pulse results: A race condition free output is shown on the left figure. The right figure is the output of an implementation with a race condition causing a noise in the results.
6.1 Development and Test Environment
6.1 Development and Test Environment
Figure 6.1: Test-case:state at t=0 left, state at t=50 right To assess the potential of the DG method a benchmark problem (Category 4, Problem 1) from the first workshop on benchmark problems in computational aero-acoustics is solved [40]. This test case consists of the propagation of acoustic waves generated by a Gaussian pulse. Since the exact solution to the Linearized Euler Equations (LEEs) is available here, the results of simulations are used to study the accuracy of the spatial discretization scheme [29].
Figure 6.1: Test-case:state at t=0 left, state at t=50 right To assess the potential of the DG method a benchmark problem (Category 4, Problem 1) from the first workshop on benchmark problems in computational aero-acoustics is solved [40]. This test case consists of the propagation of acoustic waves generated by a Gaussian pulse. Since the exact solution to the Linearized Euler Equations (LEEs) is available here, the results of simulations are used to study the accuracy of the spatial discretization scheme [29].
Table 6.2 shows the estimated error Lzerror [42] of the numerical solution computed on the GPU. The error decreases with the increase of the mesh resolution from 8xh to 1xh. The estimated error decreases also with the increase of the polynomial order. Figure 6.2 shows the error estimate Lzerror for different polynomial order p plotted against the computational time. The error estimate presents a similar trend to the reference one achieving a certain accuracy level. The original CPU implementation generates the same error estimate plot with the aberration in the third polynomial order (cf. [29]).
Table 6.2 shows the estimated error Lzerror [42] of the numerical solution computed on the GPU. The error decreases with the increase of the mesh resolution from 8xh to 1xh. The estimated error decreases also with the increase of the polynomial order. Figure 6.2 shows the error estimate Lzerror for different polynomial order p plotted against the computational time. The error estimate presents a similar trend to the reference one achieving a certain accuracy level. The original CPU implementation generates the same error estimate plot with the aberration in the third polynomial order (cf. [29]).
Table 6.4: Speedup of the eDisco on GPU and different number of CPU cores combination
Table 6.4: Speedup of the eDisco on GPU and different number of CPU cores combination
Figure 6.4: Speedup high mesh resolution Figure 6.3: Speedup medium mesh resolution
Figure 6.4: Speedup high mesh resolution Figure 6.3: Speedup medium mesh resolution
Figures 6.3 to 6.6 illustrate the reached speedups of multi-cores CPUs and the Tesla and Kepler GPUs with the one-core performance as reference. The performance measurement is done on medium and fine resolution (22300 and 99900 elements). The GPU makes a peak speedup of 66x in the medium resolution and only 50x peak performance in the fine resolution (see Figure 6.5 and 6.6).
Figures 6.3 to 6.6 illustrate the reached speedups of multi-cores CPUs and the Tesla and Kepler GPUs with the one-core performance as reference. The performance measurement is done on medium and fine resolution (22300 and 99900 elements). The GPU makes a peak speedup of 66x in the medium resolution and only 50x peak performance in the fine resolution (see Figure 6.5 and 6.6).
Figure 6.5: Speedup of a Kepler GPU with medium mesh resolution
Figure 6.5: Speedup of a Kepler GPU with medium mesh resolution
Figure 6.6: Speedup of a Kepler GPU with medium mesh resolution
Figure 6.6: Speedup of a Kepler GPU with medium mesh resolution
Table 6.5: Speedup of the eDisco on different GPUs with the 16xcores Xeon workstation as reference and threads so much as a multiprocessors register and local memory can manage. The same configuration is run on the Kepler architecture. The programmer can take more advantage from the Kepler processor with the optimization of the started number of threads and blocks and also the use of the new feature as the L2 cache.
Table 6.5: Speedup of the eDisco on different GPUs with the 16xcores Xeon workstation as reference and threads so much as a multiprocessors register and local memory can manage. The same configuration is run on the Kepler architecture. The programmer can take more advantage from the Kepler processor with the optimization of the started number of threads and blocks and also the use of the new feature as the L2 cache.
On the other hand the execution time per DOF decreases with the increase of the mesh resolution. A better resolution is coupled with a higher number of elements. The increase of the elements number let the GPU compute more elements in parallel. The GPU starts a high number of threads and if the elements number is lower than the started threads some threads are inactive, which means a computation power portion is not used, some cores are idle.
On the other hand the execution time per DOF decreases with the increase of the mesh resolution. A better resolution is coupled with a higher number of elements. The increase of the elements number let the GPU compute more elements in parallel. The GPU starts a high number of threads and if the elements number is lower than the started threads some threads are inactive, which means a computation power portion is not used, some cores are idle.
Listing A.1: Access of the Device Properties Property of the graphics processing unit device are saved in the structure cudaDeviceProp in the SDK header file (driver_types.h). The structure contains the device name, device computation power, warp size and many other device information. The listing A.1 shows the use of the API function cudaGetDeviceProperties, which accesses the device properties.
Listing A.1: Access of the Device Properties Property of the graphics processing unit device are saved in the structure cudaDeviceProp in the SDK header file (driver_types.h). The structure contains the device name, device computation power, warp size and many other device information. The listing A.1 shows the use of the API function cudaGetDeviceProperties, which accesses the device properties.
\ comparison with the error estimate of the numerical solution computed on 1xcore, xcores and 16xcores CPU and the Tesla GPU is shown in Table A.2. The error estimate f different architecture are from the same order of magnitude. The computation is done n single precision.
\ comparison with the error estimate of the numerical solution computed on 1xcore, xcores and 16xcores CPU and the Tesla GPU is shown in Table A.2. The error estimate f different architecture are from the same order of magnitude. The computation is done n single precision.
Table A.3: Speedup of the GPU and different number of CPU processors (side bump, 5900 cells)
Table A.3: Speedup of the GPU and different number of CPU processors (side bump, 5900 cells)
Table A.4: Speedup of the GPU toward different number of CPU processors(side bump, 22300 cells)
Table A.4: Speedup of the GPU toward different number of CPU processors(side bump, 22300 cells)
Related topics:
Discontinuous Galerkin MethodsEuler Equations
580 California St., Suite 400
San Francisco, CA, 94104
© 2025 Academia. All rights reserved