The Kramers-Klein equation is a special Fokker-Planck equation describing the brownian motion in a potential. This equation is used in di erent areas of biology, chemistry and physics. In this paper we present the comparison between numerical and analytical solutions in case of a linear force and usage of multi-core CPU and GPU for solving Kramers-Klein equation with di erent boundary conditions. It was shown that the GPU calculations under re ecting boundary conditions are slower as compared with other boundary conditions because in the rst case we need to call two kernels sequentially.
CUDA
Brownian motion(BM) is the random motion of a particle in a uid. Hereafter we're discussing the one-dimensional case. BM can be mathematically described in two ways. The rst one is a Langevin equation, a stochastic di erential equation for the particle's coordinate:
t) (t) >= 0 (1) here m is the mass of the particle, is the damping constant, f 0(x) is the stationary potential force per mass due to the potential f (x), (t) { deltacorrelated Gaussian random process. It's clearly seen from (1) that x(t) and v(t) = x_ (t) are stochastic quantities due to (t) in the equation. Therefore we may ask for probability density function (PDF)( (x; v; t)) of nding particle in a coordinate x with a velocity v at time t. An equation for PDF in case of BM in an external potential eld is called Kramers-Klein (KK) equation: Despite intensive analytical studies of KK equation [2], analytic solution can be found only in special cases. This makes numerical solutions as a signi cant resource for studying stochastic processes. Di culties of studying KK equation numerically come from the fact, that this partial di erential equation is of a mixed type. The problem behaves like a parabolic equation in the v-direction and like a hyperbolic equation in the x-direction. Moreover, we have di erent types of hyperbolic equations in the x-directions for v > 0 and v < 0. We used a numerical method, computational domain and the uniform grid as in [3].Central di erences scheme were used for the discretization in v - direction of rst and second derivatives, for x - direction the upwind scheme discretization was used, for t - direction forward Euler scheme was used. It is an explicit method that e ciently deals with di erent types of boundary conditions and potential elds. Convergence and stability of the method is discussed in Ref.[3]. The accuracy of the method is de ned as O( x + v2 + t). 3
Say we have a rectangular 2D computational domain(in x and v direction, xL x xR,vL v vR) with uniform grid. In the v-direction the problem doesn't have boundaries and we assume that computational domain is so large that we can assume the BC as follows: (x; vL; t) = (x; vR; t) = 0 (3) for all x; t.
In the x-direction we consider three types of BC: { Dirichlet. In analogy with v-direction we assume that computational domain is so large that particle can't be in boundaries of domain: (xL; v; t) = (xR; v; t) = 0 (4) for all v; t. { Absorbing. In this case we only have physical boundaries de ned at xL for v < 0 and at xR for v > 0: The physical meaning of the absorbing BC is that there is no particle ux from the boundary into the interior. { Re ecting. In this case we have the BC as follows: (xL; v; t) = 0; v > 0 (xR; v; t) = 0; v < 0 (xL; v; t) = (xL; v; t) (xR; v; t) = (xR; v; t) The physical meaning of the re ecting BC is that no particle would be lost at the boundary and all particles will be re ected back into the interior. (5) (6)
Parallel CPU calculations were implemented using OpenMP, the GPU calculations were done using CUDA. The numerical method we used [3] is explicit, the PDF inside the domain is fully calculated from PDF at previous time, as a consequence, calculation of PDF in each interior's node is independent of another nodes at the same time: (x; v; ti+1) = (x; v; (x; v; ti)) (7) hence the method can be easily parallelized using \stencil" parallel pattern.
In case of CPU parallelization, the boundaries were calculated separately and without parallelization because the highest performance was shown in this case. In case of GPU parallelization the boundaries were calculated in one CUDA kernel with interior in case of the Dirichlet and absorbing BC.
In case of the re ecting BC two CUDA kernels run sequentially: the rst one calclulates interior and half of boundaries in x-direction with v > 0 and the second one calculates half of the boundaries in x-direction with v < 0. We did it because in this case PDF on the boundary depends on PDF in another points at the same time(see (6)). As a result, GPU calculations with re ecting BC are slower than with the Dirichet or absorbing ones.
In CUDA calculations only global memory was used. 5
Program for CPU was written in C++ and compiled with gcc (g++ -O3) compiler. GPU implementation is also in C++ and CUDA part was compiled with proprietary nvcc compiler while g++ has been used for the host code. CPU calculations were done on PC(4x Intel Core i5 3.2 GHz, NVIDIA GTX 750), GPU calculations were done on PC and Uran supercomputer(NVIDIA Tesla M2050) from Institute of Mathematics and Mechanics UrB RAS. Calculations were done using double precision oating point numbers.
First, we compared results of the numerical solution with the analytic one in case of a linear force (parabolic potential). The analytic solution in this case was taken from [2] and is a two-variable Gaussian distribution of coordinate and velocity with the time-dependent means and variances. The error was calculated as a discrete norm:
0 x v X( (xi; vj ; t) i;j
10:5 a(xi; vj ; t))2A (8) where a is an analytic solution. In Figure 1 the numerical and analytic PDFs are shown at t = 12 (dimensionless units). Number of nodes in the x-direction is 2049 and in the v-direction is 129. It is clearly seen that numerical solution does approximate the analytic one quite accurately. 4 0.05
Then we looked on the performance of parallel CPU and GPU calculations comparing with the serial CPU ones with di erent BC. Calculations were done on an uniform grid with 2049 nodes in x-direction and with 257 nodes in vdirection, total of 526593 nodes. Time step was 1:4e 4 and calculations were stopped when t = 10. The results are shown in the Table 1. Parallel solving of KK equation provides a signi cant speed up of the calculations even in case of the parallel CPU solution. KK equation is a partial di erential equation of mixed type which describes the time evolution of a PDF of a brownian particle in an external potential. In literature the Kramers problem( nd the escape rate of a Brownian particle from a potential well) is usually solved analytically using some assumptions(harmonic approximation near minima [1]), piecewise-linear approximation [4] and high energy barriers approximation[5].
Finding PDF is a much harder task and can be solved analytically only in a few cases. In this paper we solved the KK equation numerically. Solving KK equation numerically allows us to get a solution for both nite and in nite geometry and arbitrary potential eld. We showed that both CPU and GPU parallel computing provide a signi cant gain in the solving time. Even home PC's GPU can solve the equation for an acceptable time.
BC type is important from a physical viewpoint on the problem when it has a nite geometry. Di erent types of collisions of brownian particle with a wall: the elastic(re ecting BC) and inelactic (absorbing BC) ones in uence the PDF very signi cantly[3].
From a viewpoint of computations, in case of the re ecting BC, the GPU parallelization is less e ective because in one iteration we need to call two kernels sequentially and Table 1 shows quite clearly, that the GPU acceleration is smaller in case of the re ecting BC.
From a viewpoint of computations, further works should include more accurate and e cient schemes for solving the KK equation and its parallelization[7]. Further physical researches should include numerical studies of the KK equation for the ion channel kinetics[6] with di erent physically reasonable potentials, in nite and nite geometries.