Running times per iteration of

## 1. Introduction

The efficient solution of the Frank-Kamenetskii partial differential equation through the implementation of parallelized numerical algorithms or GPUs (Graphics Processing Units) in MATLAB is a natural progression of the work which has been conducted in an area of practical import. There is an on-going interest in the mathematics describing thermal explosions due to the significance of the applications of such models - one example is the chemical processes which occur in grain silos. Solutions which pertain to the different geometries of such a physical process have different physical interpretations, however in this chapter we will consider the Frank-Kamenetskii partial differential equation within the context of the mathematical theory of combustion which according to Frank-Kamenetskii [16] deals with the combined systems of equations of chemical kinetics and of heat transfer and diffusion. A physical explanation of such a system is often a gas confined within a vessel which then reacts chemically, heating up until it either attains a steady state or explodes.

The focus of this chapter is to investigate the performance of the parallelization power of the GPU vs. the computing power of the CPU within the context of the solution of the Frank-Kamenetskii partial differential equation. GPU computing is the use of a GPU as a co-processor to accelerate CPUs (Central Processing Units) for general purpose scientific and engineering computing. The GPU accelerates applications running on the CPU by offloading some of the compute-intensive and time consuming portions of the code. The rest of the application still runs on the CPU. The reason why the application is seen to run faster is because it is using the extreme parallel processing power of the GPU to boost performance. A CPU consists of 4 to 8 CPU cores while the GPU consists of 100s of smaller cores. Together they operate to crunch through the data in the application and as such it is this massive parallel architecture which gives the GPU its high compute performance.

The methods which will be investigated in this research are implicit methods, such as the Crank-Nicolson method (

The second hurdle is the nonlinear source term which may be dealt with using different techniques. In this chapter we will implement the Newton method which acts as an updating mechanism for the nonlinear source term and in so doing maintains the implicit nature of the scheme in a consistent fashion. While the incorporation of the Newton method leads to an increase in the computation time for the Crank-Nicolson difference scheme (see [22]) there is also an increase in the accuracy and stability of the solution. As such we find that the algorithms we are attempting to employ in the solution of this partial differential equation would benefit from the processing power of a GPU.

In this chapter we will focus on the implementation of the Crank-Nicolson implicit method, employed with and without the Newton method, and two Rosenbrock methods, namely ROS2 and ROWDA3. We consider the effectiveness of running the algorithms on the GPU rather than the CPU and discuss whether these algorithms can in fact be parallelized effectively.

## 2. Model

The steady state formulation of the equation to be considered in this chapter was described by Frank-Kamenetskii [16] who later also considered the time development of such a reaction. The reaction rate depends on the temperature in a nonlinear fashion, generally given by Arrhenius’ law. This nonlinearity is an important characteristic of the combustion phenomena since without it the critical condition for inflammation would disappear causing the idea of combustion to lose its meaning [16]. Thus, in the case of a thermal explosion, the Arrhenius law is maintained by the introduction of the exponential term which acts as a source for the heat generated by the chemical reaction. As such we are able to write an equation modelling the dimensionless temperature distribution in a vessel as

where

The value of the Frank-Kamanetskii parameter [16]

where

In this chapter we consider numerical solutions for equation (1) modelling a thermal explosion within a cylindrical vessel, i.e.

where

Frank-Kamenetskii [16] obtained a steady state solution to this problem with

In Harley [21] a Crank-Nicolson- and hopscotch scheme were implemented for equation (1) subject to (4) and (5) where

In recent work (see [22]) the Crank-Nicolson method was implemented with the Newton iteration as done by Britz et al. [9] by computing a correction set in each iteration to obtain approximate values of the dependent variable at the next time step. The efficiency of the Crank-Nicolson scheme, hopscotch scheme (both of these methods were implemented with an explicit and then an implicit discretisation of the source term) and two versions of the Rosenbrock method were compared [22]. Using the pdepe function in MATLAB and the steady state solution obtained by Frank-Kamenetskii [16] as a means of comparison, it was found that the incorporation of the Newton method for the Crank-Nicolson- and hopscotch scheme led to increased running times as

Furthermore, it was shown that while the Crank-Nicolson- and hopscotch method (with or without the implementation of the Newton method) performed well in terms of accuracy for

As a consequence of these findings and due to the complexity created by the nonlinear source term which serves a critical function in the model, further work regarding faster algorithms for the solution of such an equation are of interest. This chapter will not consider the hopscotch scheme directly as an appropriate method for the solution of the Frank-Kamenetskii partial differential equation due to work done by Feldberg [15] which indicated that for large values of

## 3. Executing MATLAB on a GPU

The advantage of using the Parallel Computing Toolbox in MATLAB is the fact that it allows one to solve computationally and data-intensive problems using multicore processors, GPUs, and computer clusters. In this manner one can parallelize numerical algorithms, and in so doing MATLAB applications, without CUDA or MPI programming. Parallelized algorithms such as parfor, used within the context of what is usually a for loop, allows you to offload work from one MATLAB session (the client) to other MATLAB sessions, called workers. You can use multiple workers to take advantage of parallel processing and in this way improve the performance of such loop execution by allowing several MATLAB workers to execute individual loop iterations simultaneously. In this context however we are not able to implement in-built MATLAB functions such as parfor due to the numerical algorithms which we have chosen to consider. The

The Parallel Computing Toolbox in MATLAB allows one to create data on and transfer it to the GPU so that the resulting GPU array can then be used as an input to enhance built-in functions that support them. The first thing to consider when implementing computations on the GPU is * keeping* the data on the GPU so that we do not have to transfer it back and forth for each operation - this can be done through the use of the gpuArray command. In this manner computations with such input arguments run on the GPU because the input arguments are already in the GPU memory. One then retrieves the results from the GPU to the MATLAB workspace via the gather command. Having to recall the results from the GPU is costly in terms of computing time and can in certain instances make the implementation of an algorithm on the GPU less efficient than one would expect. Furthermore, the manner in which one codes algorithms for GPUs is of vital importance given certain limitations to the manner in which functions of the Toolbox may be implemented (see [25]). More importantly however, is whether the method employed can allow for the necessary adjustments in order to improve its performance. In this chapter we will see that there are some problems with implementing the kind of algorithms considered here on the GPU.

In this chapter we are employing MATLAB under Windows 7 (64 bits) on a PC equipped with an i7 2.2 GHz processor with 32 GB of RAM.

### 3.1. Crank-Nicolson implicit scheme

We will implement the Crank-Nicolson method while maintaing the explicit nature of the nonlinear source term and also apply the method by computing a correction set in each iteration to obtain approximate values of the dependent variable at the next time step through the use of the Newton method [26]. The methodology will be explained briefly here; the reader is referred to [7, 8, 9] for clarification.

When implementing the Crank-Nicolson method we employ the following central-difference approximations for the second-and first-order spatial derivatives respectively

while a forward-difference approximation

is used for the time derivative. We implement a Crank-Nicolson scheme by approximating the second-derivative on the right-hand side of (1) by the implicit Crank-Nicolson [12] approximation

In a similar fashion the first-derivative on the right-hand side becomes

To impose zero-shear boundary conditions at the edges we approximate the spatial first-derivative by the central-difference approximation (7) which leads to the following condition

As mentioned before the boundary condition (5a) at

which simplifies the equation for the case

to equation (1) at

as the initial difference scheme with

where

As indicated by the boundary conditions (5a) and (5b) we consider the problem for

and is solved as follows

The inverse of

### 3.2. Rosenbrock method

We now consider two particular Rosenbrock methods, ROS2 and ROWDA3, as a means of comparison for the effectiveness of the methods discussed in the previous section. The Rosenbrock methods belong to the class of linearly implicit Runge - Kutta methods [11, 17]. They were used successfully for the numerical solution of non-electrochemical stiff partial differential equations, including equations of interest to electrochemistry. For further information regarding the particulars of such methods interested readers are referred to the numerical literature of [17, 28, 29, 30].

The reason for the use of the Rosenbrock methods in this paper is the ease with which they are able to deal with the nonlinear source term and the fact that no Newton iterations are necessary. The advantages of these methods are great efficiency, stability and a smooth error response if ROS2 or ROWDA3 are used (see [4] for instance) and the ease with which they are able to handle time-dependent and/or nonlinear systems.

We consider equation (1) as the following system

which along with

where

such that

In order to implement the Rosenbrock methods a number

where we define

where the

In this chapter we implement the ROS2 and ROWDA3 methods though there are other variants of the Rosenbrock methods. Lang [24] described a L-stable second-order formula called ROS2. A third-order variant thereof is called ROWDA3 and described by Roche [28] and later made more efficient by Lang [23]. The latter is a method favoured by Bieniasz who introduced Rosenbrock methods to electrochemical digital simulation [4, 5]. For a more detailed explanation and discussion regarding the method and its advantages refer to [7]. The focus of the work done here is with regards to whether the Rosenbrock algorithms lend themselves toward parallelized implementation. It has already been noted that functions such as the parfor command cannot be used in this instance. It now remains to consider the method’s performance when run on a GPU via the MATLAB Parallel Computing Toolbox.

## 4. Discussion of numerical results

The results noticed, as per Table 1, indicate the extent to which implementing the code on the GPU slows down overall performance of the

To gain a better understanding of these results we consider the baseline structure for our

In doing so, we realise that the main components thereof are matrix and elementwise vector operations. In order to understand why we are achieving the results we do (see Table 1) we run a few simple ’test’-codes to consider the speed taken by the CPU vs. the GPU to perform such elementary operations as

The next step in this evaluative process is to now consider the speed at which vector operations are performed. This was done in a similar fashion to the previous case by considering the speed taken by the CPU and GPU to perform the elementwise operation

CPU | 7.8449e-05 | 2.7737e-04 | 2.3002e-05 | 9.1288e-05 |

GPU | 0.0014 | 0.0960 | 0.0033 | 0.0063 |

CPU | 1.4783e-05 | 2.0018e-04 | 1.5354e-05 | 6.3574e-05 |

GPU | 8.0940e-04 | 0.0047 | 0.0028 | 0.0047 |

CPU | 2.4573e-05 | 0.0033 | 1.7146e-05 | 6.8119e-05 |

GPU | 0.0048 | 0.4731 | 0.0042 | 0.0073 |

.0137 | 0.0025 | 0.0059 | 0.0138 |

.0133 | 0.0024 | 0.0067 | 0.0149 |

.0146 | 0.0025 | 0.0057 | 0.0128 |

.0145 | 0.0024 | 0.0059 | 0.0133 |

.0146 | 0.0012 | 0.0064 | 0.0150 |

.0275 | 0.0026 | 0.0061 | 0.0160 |

.0151 | 0.0030 | 0.0062 | 0.0151 |

.0160 | 0.0042 | 0.0063 | 0.0140 |

As a means of further investigation we consider the

### 4.1. Influence of changing parameter values on the running time of the algorithms

Just a few brief comments upon the results obtained for

## 5. Concluding remarks

The implementation of numerical algorithms such as those considered in this chapter are widely used for the solution of many differential equations which model physical processes and applications. As such it is of vital importance that we be able to perform such calculations at high speed given the requirement of fine grids to improve accuracy. It is in this context that the use of GPUs becomes of prime importance. However it is not simply a matter of running the algorithm on the GPU - the method employed needs to lend itself to being adjusted in the required manner so that the parallel processing power of the GPU may be taken advantage of. Though we found that the numerical methods considered here were not entirely suited to being implemented on the GPU as we would have hoped we were able to explain why this was the case.

This work has investigated the effectiveness of matrix and elementwise operations when run on a GPU vs. a CPU and found that the speed taken to do such operations heavily relies on the choice of

We also discovered the power of the in-built function arrayfun which was able to improve upon the performance of the GPU with regards to computing time to the degree that it outperformed the CPU even for a grid with ’large’

What the work in this chapter has shown is that the structures of the GPU and the Parallel Computing Toolbox in MATLAB are such that while certain algorithms have the ability to be adjusted for improved performance not all methods do. In particular it seems clear that implicit methods with matrix and vector operations will in fact run much slower on the GPU than the CPU. Thus whether GPU computing is able to improve the performance of a numerical scheme is very much dependent upon the type of computations which need to be done. In our case we discovered that the implicit and nonlinear nature of our numerical schemes do not lend themselves towards improved performance via the implementation of the parallel processing power of a GPU.

### Acknowledgement

I would like to thank Mr. Dario Fanucchi for invaluable discussions.

## References

- 1.
Abd-Salam El M. R. Shehata M. H. 2005 The numerical solution for reaction diffusion combustion with fuel consumption Appl. Math. Comp., 160:423U˝-435. - 2.
Anderson C. A. Zienkiewicz O. C. 1974 Spontaneous ignition: finite element solutions for steady and transient conditions, J. Heat Transfer,96 3 398 404 - 3.
Anderson D. Hamn´ en. H. Lisak M. Elevant T. Persson H. 1991 Transition to thermonuclear burn in fusion plasmas Plasma Physics and Controlled Fusion,33 10 1145 1159 - 4.
Bieniasz L. K. 1999 Finite-difference electrochemical kinetic simulations using the Rosenbrock time integration scheme - 5.
Bieniasz L. K. Britz D. 2001 Chronopotentiometry at a Microband Electrode: Simulation Study Using a Rosenbrock Time Integration Scheme for Differential-Algebraic Equations, and a Direct Sparse Solve r, - 6.
Boddington T. Gray P. 1970 Temperature profiles in endothermic and exothermic reactions and interpretation of experimental rate data, Pr oc. Roy. Soc. Lond Ser A-Mat. Phys. Sci.,320 1540 71 100 - 7.
Britz D. 2005 Digital Simulation in Electrochemistry, 3rd Edition, Lecture Notes in Physics, Springer, 3 − 540 −23979 − 0, Berlin Heidelberg - 8.
Britz D. Baronas R. Gaidamauskaitė E. Ivanauskas F. 2009 Further Comparisons of Finite Difference Schemes for Computational Modelling of Biosensors, Nonlinear Analysis: Modelling and Control,14 4 419 433 - 9.
Britz D. Strutwolf J. &Østerby O. 2011 Digital simulation of thermal reactions, Appl. Math. and Comp., 218(4), 15:1280-1290 - 10.
Chambré P. L. 1952 On the solution of the Poisson-Boltzmann equation with application to the theory of thermal explosions, J. Chem. Phys.,20 1795 1797 - 11.
Chan Y. N. I. Birnbaum I. Lapidus L. 1978 Solution of Stiff Differential Equations and the Use of Imbedding Techniques Ind. Eng. Chem. Fundam.,17 3 133 148 - 12.
Crank J. Nicolson E. 1947 A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type Proc. Camb. Phil. Soc.,43 50 67 - 13.
Crank J. Furzeland R. M. 1977 The treatment of boundary singularities in axially symmetric problems containing discs J. Inst. Math. Appl.,20 3 355 370 - 14.
Evans D. J. Danaee A. 1982 A new group Hopscotch method for the numerical solution of partial differential equations SIAM J. Numer. Anal.,19 3 588 598 - 15.
Feldberg S. W. 1987 Propagational inadequacy of the hopscotch finite difference algorithm: the enhancement of performance when used with an exponentially expanding grid for simulation of electrochemical diffusion problems J. Electroanal. Chem.,222 101 106 - 16.
Frank-Kamenetskii D. A. 1969 Diffusion and Heat Transfer in Chemical Kinetics Plenum Press, New York - 17.
Hairer E. Wanner G. 1991 Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems, Springer-Verlag, 3 − 540 − 60452 − 9, Berlin - 18.
Harley C. Momoniat E. 2007 Steady state solutions for a thermal explosion in a cylindrical vessel - 19.
Harley C. Momoniat E. 2008 Instability of invariant boundary conditions of a generalized Lane-Emden equation of the second-kind - 20.
Harley C. Momoniat E. 2008 Alternate derivation of the critical value of the Frank-Kamenetskii parameter in the cylindrical geom etry, - 21.
Harley C. 2010 Explicit-implicit Hopscotch method: The numerical solution of the Frank-Kamenetskii partial differential equation, Journal of Applied Mathematics and Computation,217 8 4065 4075 - 22.
Harley C. 2011 Crank-Nicolson and Hopscotchmethod: An emphasis onmaintaining the implicit discretisation of the source term as a means of investigating critical parameters. Special Issue on’Nonlinear Problems: Analytical and Computational Approach with Applications’, - 23.
Lang J. 1996 High-resolution self-adaptive computations on chemical reaction-diffusion problems with internal boundaries - 24.
(Lang J. 2001 ). Adaptive Multilevel Solution of Nonlinear Parabolic PDE Systems, Springer, 9783540679004, Berlin - 25.
1994 201 The MathWorks, Inc. ©1994-2012. Parallel Computing Toolbox Perform parallel computations on multicore computers, GPUs, and computer clusters, http : //www.mathworks.com/products/parallel − computing - 26.
Press W. H. Teukolsky S. A. Vetterling W. T. Flannery B. P. 1986 Numerical Recipes in Fortran, 2nd Edition, Cambridge University Press, 0 − 521 − 43064 − X, Cambridge - 27.
Rice O. K. 1940 The role of heat conduction in thermal gaseous explosions, J. Chem. Phys.,8 9 727 733 - 28.
Roche M. 1988 Rosenbrock methods for differential algebraic equations Numerische Mathematik,52 45 63 - 29.
Rosenbrock H. H. 1963 Some general implicit processes for the numerical solution of differential equations - 30.
Sandu A. Verwer J. G. Blom J. G. Spee E. J. Carmichael G. R. 1997 Benchmarking Stiff ODE Solvers for Atmospheric Chemistry Problems II: Rosenbrock Solvers, Atmospheric Environment,31 3459 3472 - 31.
Steggerda J. J. 1965 Thermal stability: an extension of Frank-Kamenetskii’s theory, J. Chem. Phys.,43 4446 4448 - 32.
Zhang G. Merkin J. H. Scott S. K. 1991 Reaction-diffusion model for combustion with fuel consumption: Ii. Robin boundary conditions IMA J. Appl. Math.,51 69 93 - 33.
Zhang G. Merkin J. H. Scott S. K. 1991 Reaction-diffusion model for combustion with fuel consumption: I. Dirichlet boundary conditions IMA J. Appl. Math.,47 33 60 - 34.
Zeldovich Y. B. Barenblatt G. I. Librovich V. B. Makhviladze G. M. 1985 The Mathematical Theory of Combustion and Explosions Consultants Bureau, New York