Implementation of Numerical Methods for Partial Differential Equation Using Parallel Computing

ISSN (Online) 2394-6849 International Journal of Engineering Research in Electronics and Communication Engineering (IJERECE) Vol 4, Issue 8, August 2017 Implementation of Numerical Methods for Partial Differential Equation Using Parallel Computing [1] [3] [4] Mr. Akshay A. Jadhav, [2] Mrs. Trupti P. Agarkar, Mr. Vishwesh A. Vyawahare, Mr. Mukesh D. Patil [1][2][3] Department of Electronics Engineering, [4]Department of Electronics and Telecommunication Engineering [1][2][3][4] Ramrao Adik Institute of Technology Abstract - The performance and use of parallel computing in the field of differential calculus is increased tremendously opening up new avenues for applying these in the field of numerical computation for high speed performance. The computation time required to find analytical as well as numerical solution is tested and compared. In this work we have harnessed this property of GPU to accelerate the grid point calculations for numerical calculations and the performance of numerical method using CPU and GPU is compared. The numerical Methods for integer order PDE are studied, analyzed and implemented on GPU using parallel computing toolbox of MATLAB. The finite difference methods of PDE like explicit, implicit method are tested for the results, for parabolic, hyperbolic and elliptical type of PDE’s. The positive speed up is achieved for elliptical type of PDE. The verification of results with the analytical solution is made by the mean square error. Index Terms: Parallel Computing Toolbox, GPU Computing, Partial Differential Equation, Numerical Method. I.INTRODUCTION which typically handles computation only for graphics. It used to perform computation in the Partial differential equations (PDEs) are used to applications traditionally handled by the central describe most of the physical phenomena‟s in the processing unit (CPU). The use of multiple video domains like fluid dynamics, electricity, magnetism, cards in one computer, or large numbers of graphics mechanics, optics or heat flow. To solve these PDE‟s chips, further parallelizes the already parallel nature of it‟s difficult to do the pen and paper calculations or to graphics processing. Addition to this even a single use the analytical methods. Finding the numerical GPUCPU framework provides advantages that solution or the approximate solution over the multiple CPUs on their own can‟t offer due to analytical solution is preferably considered. The specialization in each chip [3]. To solve the hyperbolic numerical solutions are extremely abundant, because (Wave Equation) and elliptical (Poisson‟s Equation) sometimes we don‟t have an analytical approach or the type of PDE we used the explicit method as the finite analytical solution is too slow. Instead of computing difference numerical method. for 15 hours and getting an exact solution, we rather compute for 15 seconds and get a good approximation. II. INTRODUCTION TO PDE While solving the partial differential equations for PDE (Partial Differential Equation) is a partial example regularly used diffusion equation Ut = derivative of he equation and contains an unknown Uxx(1D) or Ut = Uxx + Uyy(2D) we require a rigorous function depending on more than one variable. For a calculations, to overcome this we use a GPGPU general 2nd order linear PDE AUxx+ BUxy +CUyy pipeline. It is a process of analyzing a data in terms of +DUx + EUy+FU = G [1], [2] Where A, B, C, D, E, F, graphical arrays and does the parallel processing in G are all functions of x,y. This equation is classified one or more GPU and CPU. GPUs are having large w.r.t sign of the discriminate ∆S i.e. ∆S = B2 ∆ 4AC. number of cores but at the same time they operate at The few regularly tested PDE‟s are classified based on lower frequencies. Thus, the speed of operation of the values of ∆ [1]. GPU in comparison with CPU is greater for pictures and graphical data. Transferring data into graphical The PDE is solved to find the unknown function U. form and then with the help of GPU, scan and analyze This solution for PDE can be of two types first is the data will result in large speedup. General-purpose analytical (i.e. Exact) solution and second is the computing on graphics processing units (GPGPU) numerical (i.e. Approximate) solution, Which will means the use of a graphics processing unit (GPU), satisfies any initial and boundary conditions mentioned with PDE [7]. Analytical solutions are the All Rights Reserved © 2017 IJERECE 29  ISSN (Online) 2394-6849 International Journal of Engineering Research in Electronics and Communication Engineering (IJERECE) Vol 4, Issue 8, August 2017 exact solutions and can be achieved by pen and paper calculations where as Numerical solutions are approximate solutions and we are unable to find it in finite time and it is difficult to solve it by using pen and paper. This distinction between analytical and numerical can vary with the applications. In Most of the applications PDEs do not have analytical solutions so it is a necessary to use a numerical procedure to find an approximate solution [3]. To find an approximate solution we have to approximate the independent variables at discrete values and implement the numerical method via a computer program. To find the unique solution we need to concentrate on three things 1) Domain Dependency: The domain for PDE means Fig. 1. Initial and Boundary Condition considered the region in which PDE is applied in space and time for grid Solution [1] variables. In practical applications, PDE will be defined on a bounded domain, or in a closed boundary III. PARALLEL COMPUTING TO SOLVE [4], [6]. NUMERICAL METHODS FOR PDE 2) Initial Condition: we have to define a initial While dealing with big technical problem we came condition such that it specifies the value of unknown across the issues like long computational running time function u at the initial time on the domain [1], [2], for the process or the part of the code is [11]. computationally intensive or the system of application for which we are coding has large input data set. All 3) Boundary Condition: these issues hammer the speed up of process and slow Dirichlet Boundary Condition: The values of down the performance. The best solution to deal with unknown function U on the domain boundary are this problem is parallel processing. specified by the dirichlet boundary condition. General purpose computation on GPUs (GPGPU) uses Numerically, this condition guarantees the existence a graphics processing unit (GPU) for general-purpose and uniqueness of the PDE [1], [11]. parallel processing applications. No dependencies or communications between tasks is the most important Von-Neumann Boundary condition: The von- criteria for the parallel computing [14].Vectorization is Neumann condition is defined for the domain one of the concern thing for GPU coding, Most boundary as the normal derivatives with respect to a operations on the GPU operate in a vectorized fashion, variable. The von-Neumann condition can be defined one operation can be performed on up to four values at by the condition: (y) [2], [6], [11]. once. This functionality of producing output from two Unlike the Dirichlet condition, the Neumann different independent input data sources is useful in condition ensures only that the equation has a unique graphics because almost every basic data type is a solution up to an additional constant. When the vector (either 2-D, 3-D, or 4-Dimensional). [12], [15] Neumann condition is constant, it implies that the rate When MATLAB‟s GPU enabled functions are of variation of u across the boundary does not change executed on the GPU, data will be transferred from [1], [2]. MATLAB‟s workspace to GPU device memory. The command „gpuArray‟ provides a specific array type for such data transfer, then GPU enabled functions can run on these data. The command „gather‟ returns the calculated results, which are stored in GPU, back to MATLAB‟s workspace [14].The procedure of GPU computing in MATLAB is as shown in figure(2) All Rights Reserved © 2017 IJERECE 30  ISSN (Online) 2394-6849 International Journal of Engineering Research in Electronics and Communication Engineering (IJERECE) Vol 4, Issue 8, August 2017 process [14]. We are using MATLAB as our software platform which supports the NVIDIA GPU, with the help of that we came to understand which part of the application code is consuming more time so, that part can be parallelized. This part is then transferred to the Fig. 2. Steps followed in GPU computing GPU for accelerated parallel processing, while the rest of the program is executed sequentially by the CPU. The advantage of MATLAB GPU programming is that The executed output data set is gathered by the GPU users can easily utilize GPU computing by adding just memory and transferred back to the CPU memory, few more commands to their original MATLAB which will be available in MATLAB workspace and codes. Disadvantages include that only limited for further calculations. As a result to all this the MATLAB functions are GPU enabled and the computationally complicated analysis of mathematical computing efficiency of MATLAB GPU is less than equations and methods programs runs faster and this those codes written in CUDA [14]. results in efficient processing for this computationally intensive codes. [12], [15] IV. GPU COMPUTING WITH MATLAB GPU accelerated computing off loads compute- intensive portions of the application to the GPU, while the remainder of the code still runs on the CPU. A simple way to understand the difference between a GPU and a CPU is to compare how they process tasks [8]. A CPU consists of a few cores optimized for sequential serial processing while a GPU has a massively parallel architecture consisting of thousands of smaller, more efficient cores designed for handling multiple tasks simultaneously [8], [14]. Fig. 3. Working of GPU Acceleration [14] we are using the hybrid GPU-CPU system in which Fig. 4. Workflow of GPU Computing CPU will be treated as a host or a primary processor and the GPU is called the device or a secondary V. PARALLEL COMPUTING TOOLBOX (PCT) processor. When we run an application code on the machine the CPU as main processor decides which is The Parallel Computing Toolbox (PCT) is a part of the “computationally-intensive” part of the program MATLAB toolbox. It helps us to solve [9]. This part actually takes large amount of computationally intensive and data-intensive problems computational time and due to which its data set can using MATLAB speedily on our local multi core be executed in parallel on GPU with the ongoing computer or on Shared Computing Cluster. Parallel processing operations such as gpuArray and gather, All Rights Reserved © 2017 IJERECE 31  ISSN (Online) 2394-6849 International Journal of Engineering Research in Electronics and Communication Engineering (IJERECE) Vol 4, Issue 8, August 2017 implements the tasks and perform data-parallel algorithms in MATLAB. Converting serial MATLAB applications to parallel MATLAB applications it generally requires few code modifications and no programming in a low-level language is necessary. We can run our parallel applications interactively or in batch [14]. Considering the Dirichlet boundary conditions and the initial conditions, and considering the mesh which is The main advantages of parallel computing are: uniformly spaced [3], [8]. 1) Significant amount of time and money are saved. 2) Larger problems can be solved. 3) Achievement of parallelism. To find the solution by explicit method for wave 4) Usage of non-local resources. equation we have to discritize it over a domain of boundaries,for that will use finite difference Parallel Computing Toolbox helps us to solve approximations to replace the second order computationally and data-intensive problems using derivatives,namely multi core processors, GPUs, and computer clusters. High-level constructs parallel for-loops, special array types, and parallelized numerical algorithms which helps to parallelize MATLAB applications without CUDA or MPI programming. The toolbox uses the full processing power of multi core desktop Without Since the errors are of orders of ( t)2 and ( x)2, we changing the code, we can run the same applications expect to be able to choose the space and time step on a computer cluster or a grid computing service. sizes of comparable magnitude: t x VI. HARDWARE ASPECTS: leopard - GPU compute node has  Tesla K40 GPU with 2880 cores 12 GB RAM.  Intel(R) Core(TM) i7-7500U CPU @ 2.90GHz with 6 cores 16 GB RAM. VII. NUMERICAL METHOD TO SOLVE HYPERBOLIC PDE The time dependent (transient) phenomena are described by Hyperbolic Equations for e.g. sea Fig. 5. schematic diagram of a five point formula for tides,String vibrations,wave motion,liquid flow. From Explicit scheme [3] the general 2nd order linear PDE for the value of Substituting the finite difference formulae into the discriminant ∆ = B2-4AC > 0 particular PDE is partial differential equation(1), and rearranging terms, classified as hyperbolic PDE. The 1D wave equation, we are led to the iterative system utt = c2uxx (c = constant > 0), is the simplest model of a physical system involving a hyperbolic PDE. [8] Let us consider the second order hyperbolic equation in for i = 1,2,3,..... and j = 1,....Nt - 1. one space dimension for the numerical approximations to the solution values at the mesh points.The positive parameter All Rights Reserved © 2017 IJERECE 32  ISSN (Online) 2394-6849 International Journal of Engineering Research in Electronics and Communication Engineering (IJERECE) Vol 4, Issue 8, August 2017 Stability Analysis: For the stability of 2, denoting and using the von Neumann technique, we get the condition as under which this method is stable [6]. Example: The code solves the above problem with string of length l, T = 1 by the Explicit scheme. The boundary and initial conditions are chosen properly such that our Fig. 7. Numerical solution for hyperbolic equation by problem has the analytic solution Explicit Scheme TABLE I RESULTS OBTAINED FOR WAVE EQUATION BY EXPLICIT SCHEME Fig. 8. Error Analysis for hyperbolic equation by Explicit Scheme 1) The Analytical Solution (figure 6) and Numerical Solution (figure 7) are almost matching with the error specified in (figure 8) result table for each space and time step combination. 2) As step size for space or time decreases that is as we are going to increase space and time samples the truncation error tends to zero and solution becomes more stable also it approaches towards the exact solution.The mean square error(MSE) is of the order Fig. 6. Analytical solution for hyperbolic equation by of 10-29 Explicit Scheme 3) Speed Up for GPU is negative the reason for this is the for loop which we are unable to vectorize contains the dependency issues. As the method itself has the dependencies on neighbouring values unless they are All Rights Reserved © 2017 IJERECE 33  ISSN (Online) 2394-6849 International Journal of Engineering Research in Electronics and Communication Engineering (IJERECE) Vol 4, Issue 8, August 2017 calculated it is difficult to find further values, but we with a homogeneous Dirichlet boundary condition are trying to overcome this. over the domain Ω = (0; 1)2. The right-hand side function f is chosen as f(x, y) = 2𝜋2sin(𝜋 x)sin(𝜋 y) 4) one important thing is both the time and space such that the exact solution is given by u(x; y) = samples are same to achieve the stability criteria. sin(πx)sin(πy). Following results were obtained on specified VIII. NUMERICAL METHOD TO SOLVE hardware: ELLIPTICAL PDE Eliptical partial differential equatins arise usually from equilibrium or steady state problems and their solutions. From the general 2nd order linear PDE for the value of discriminant ∆= B2 - 4AC < 0 particular PDE is classified as elliptical PDE. The very well known elliptical equations are poisson‟s equation, which is often written as ∇2∅=f(x,y) and TABLE II RESULTS OBTAINED FOR POISSON’S Laplace‟s equation . EQUATION BY EXPLICIT SCHEME Suppose that Ω is a bounded domain of R2 with boundary∂Ω The equation is said to be elliptic if b2 � ac < 0 for all points (x; y) Considering the Poisson‟s equation [9]. As before, we assume that is covered by a uniform grid The approximate solution at point (xi, yj) is denoted Fig. 9. Analytical solution for Eliptical equation by uji . Using the secondorder central difference [1], to approximate the derivatives uxx and uyy in 4 respectively, we obtain the five-point difference scheme: Fig. 10. Numerical solution for Eliptical equation Example: All Rights Reserved © 2017 IJERECE 34  ISSN (Online) 2394-6849 International Journal of Engineering Research in Electronics and Communication Engineering (IJERECE) Vol 4, Issue 8, August 2017 GPU depends on its number of cores and compute capability. The disadvantage for this parallel computing with MATLAB is, we don‟t have any control over the number of cores of GPU to be used for the execution. The usage of cores of GPU will be decided by the MATLAB. We will have that advantage in CUDA. To achieve speed up with the GPU our application code must satisfy the important criteria as the fact that sending the data between CPU and GPU must take less time than the performance gained by running on GPU, then only it will be good candidate of GPU functionality.The promising power of GPU with MATLAB is the main advantage. GPU enabled functions and gpuArray are useful for speedup without the use of CUDA programming. Fig. 11. Error Analysis for Eliptical equation REFERENCES Result Analysis: [1] Farlow, Stanley J. Partial differential equations for 1) The Analytical Solution (figure 9) and Numerical scientists and engineers. Courier Corporation, 1993. Solution (figure 10) are almost matching with the error (figure 11) specified in result table II for each space [2] Smith, Gordon D. Numerical solution of partial and time step combination. differential equations: finite difference methods. Oxford university press, 1985. 2) As we are going to increase space and time samples [3] Lindfield, George, and John Penny. Numerical the truncation error tends to zero and solution becomes methods: using MATLAB. Academic Press, 2012. more stable also it approaches towards the exact solution. [4] Ames, William F. Numerical methods for partial differential equations. Academic press, 2014. 3) In case of speed ups for large number of input samples better speed up (almost 3 times) is achieved [5] Peaceman, Donald W., and Henry H. Rachford, Jr. between vectorized code and sequential code. “The numerical solution of parabolic and elliptic differential equations.” Journal of the Society for 4) for the last sample we can see that at J=140 GPU Industrial and Applied Mathematics 3.1 (1955): 28-41. reaches its memory limit and through a error message as „out of memory‟.The error is Not a Number (NAN). [6] Arnold, Douglas N. “Stability, consistency, and convergence of numerical discretizations.” IX. CONCLUSION Encyclopedia of Applied and Computational Mathematics. Springer Berlin Heidelberg, 2015. 1358- To evaluate the benefits of use of the GPU to solve 1364. Partial Differential Equations we ran a benchmark study between analytical solution and the numerical [7] John, Volker. “Numerical methods for partial solution, in which we measured the amount of time the differential equations.” Lecture notes (2013). algorithm took to execute various time steps for increasing number of grid size. It is observed that [8] Trangenstein, John. “Numerical solution of increase in complexity of the functions increases the hyperbolic conservation laws.”Cambridge University speed-up. The improved speed up was achieved for Press, (2009). large number of samples. [13] [14] On the hardware side the execution speed offered by [9] Urroz, G. E. “Numerical Solution of Laplace the Equation.” Utah State University paper (2004). All Rights Reserved © 2017 IJERECE 35  ISSN (Online) 2394-6849 International Journal of Engineering Research in Electronics and Communication Engineering (IJERECE) Vol 4, Issue 8, August 2017 [10] Peter Knabner, Lutz Angermann. “Numerical Methods for Elliptic and Parabolic Partial Differential Equations” Springer-Verlag,New York Berlin Heidelberg,2000. [11] Arnold, Douglas N. “Lecture notes on Numerical Analysis of Partial Differential Equations.” (2012). [12] Chun-Hasiang Han, “Introduction to Matlab GPU Acceleration for COmputational Finance.” Natational Tsing-Hua University,Taiwan, December 26, 2013. [13] “Partial Differential Equation”. En.Wikipedia.Org, 2017, https:// en.wikipedia .org/ wiki/ Partial-differential-equation. [14] “NVIDIA On GPU Computing And The Difference Between Gpus And Cpus”. Nvidia.Com, 2017, http://www.nvidia.com/object/what-isgpu- computing.html. [15] “Graphics Processing Unit”. En.Wikipedia.Org, 2017, https:// en .wikipedia .org /wiki / Graphics- processing -unit. All Rights Reserved © 2017 IJERECE 36