We describe a parallel implementation of a di erence scheme for the advection equation with time delay on a hybrid architecture computation system. The di erence scheme has the second order in space and the rst order in time and is unconditionally stable. Performance of a sequential algorithm and several parallel implementations with the MPI technology in the C++ language has been studied.
First order partial di erential equations with time delay, also known as advection
equations, with distributed parameters arise in the modeling of birds migration,
expansion of in uence virus and transfer of nuclear particles [
The qualitative theory of partial functional di erential equations (PFDE) is
developed quite well (see, for example, [
Especially e ective di erence schemes for PFDE of parabolic, hyperbolic and
advection type were elaborated in [6{10]. The main idea in these works is a
separation principle the essence of which is the separation of nite and in nite
dimensional components in the structure of PFDE. To take into account the time
delay e ect, interpolation and extrapolation of discrete prehistory is used. This
extrapolation is needed for the realization of implicit methods and allows the
authors to avoid the necessity of solving nonlinear systems. Our approach is close
to [
A dramatic improvement of supercomputers architecture and their
performance led to the increased interest in the elaboration of parallel numerical
methods. Below we survey some of approaches. The rst class of methods of concern
are a domain decomposition (DD) methods which are based on a physical
decomposition of a global solution domain. The very rst of them is a classical
alternating Schwarz method for the elliptic equation which is sequential in its
nature. In order to obtain parallel analogs of Schwarz alternating method many
new techniques have been introduced, such as additive Schwarz methods,
parallel multilevel precondition algorithms, parallel weighted Schwarz algorithms,
parallel subspace correction methods and etc., see e.g. [
A parallel implementation of a method of the semi-Lagrangian type for the
advection equation was considered in [
Time domain decomposition algorithm for the parallel-in-time approximation
of solution of advection equation was developed in [
In this article we present a parallel version of one algorithm which is based
on a di erence scheme, rst published in [
Let > 0 and consider the following advection equation with aftere ect (1) (2) (3) + a where x 2 [x0; X] is a state and t 2 [t0; ] is time; u(x; t) is an unknown function; ut(x; ) = fu(x; t + ); < 0g is a prehistory-function of the unknown function to the moment t and a > 0 is a constant. Together with the advection equation we have the following initial and the boundary conditions u(x; t) = '(x; t); x 2 [x0; X]; t 2 [t0
; t0]; u(x0; t) = g(t); t 2 [t0; ]:
We adopt the compatibility condition g(t0) = '(x0; t0). Questions of the
existence and uniqueness of a solution to the stated boundary value problem
(1){(3) were considered in [
We denote by Q = Q[ ; 0) the set of functions u( ) that are piecewise continuous on [ ; 0) with a nite number of points of discontinuity of the rst kind and right continuous at the points of discontinuity. We de ne a norm on Q by kukQ = sup 2[ ;0) ju( )j. We additionally assume that the functional f (x; t; u; v( )) is given on [0; X] [t0; ] arguments: R
Q and is Lipschitz in the last two 9 Lf 2 R 8 x 2 [x0; X]; t 2 [t0; ]; u1 2 R; u2 2 R; v1 2 Q; v2 2 Q : f (x; t; u1; v1( )) f (x; t; u2; v2( ))
Lf ju1 u2j + kv1( ) v2( )kQ : 2
We consider an equidistant partition of [x0; X] into parts with step h = (X x0)=N and split the segment of variation of the time variable [t0; ] into parts with the step . The uniform grid fxi; tj gjM=0iN=0 can be constructed, here xi = x0 + ih; i = 0; : : : ; N; and tj = t0 + j ; j = 0; : : : ; M: Denote by uij approximations of the functions u(xi; tj ); i = 0; : : : ; N; j = 0; : : : ; M; at the nodes. Without loss of generality and to simplify the narration we assume that the value = = m is a natural number.
Since functional f (xi; tj ; u(xi; tj ); utj (xi; )) may depend on values of the function u between grid nodes, interpolation may be needed. For every xed node (xi; tj ) and time delay 2 [ ; 0) there are only two possibilities: if tj + t0; interpolation is not needed, we use the initial function, u(x(i); tj + ) = '(x(i); tj + ), otherwise we use the interpolation as described below.
For every xed time moment tj ; j = 1; : : : ; M; we introduce its discrete prehitory uik j =
i uk j maxf0; j mg k j : De nition 1. A mapping I de ned on the set of all admissible discrete prehistories and acting by the rule I : uik j ! vi;j ( ) 2 Q[ minf ; tj g; 0] is called an interpolation operator for the discrete history.
Let us give an example of a concrete interpolation operator, which has the properties required for the numerical method that we are going to construct. For the discrete prehistory uik j we de ne
1 vi;j (tj + ) = (tl+1 tj ) uli + (tj + tl) uli+1 ; tl tj + tl+1: (4) We say an interpolation operator has order of error p on the exact solution, if there exist constants C1 and C2 such that, for all i = 1; : : : ; N; j = 1; : : : ; M; and t 2 [maxf0; tj g; tj ] the following inequality holds: For example, the operator of interpolation (4) is of the second order.
We consider the following family of di erence schemes (parametrized by 0 s 1), with j = 0; : : : ; M 1: 1 uj+1
i uj+1 Because of (1) we have 1 a fj0 g_j : ; and due to (3) we obtain Theorem 1. Let the exact solution u(x; t) of problem (1) be thrice continuously di erentiable with respect to state variable x, twice continuously di erentiable with respect to time t and the rst derivative of the solution with respect to x is continuously di erentiable in t. Then if 2s > 1 method (5) converges with order h2 + ; i. e. there exists a constant C such that ku(xi; tj ) uij k C(h2 + ) for all i = 0; : : : ; N and j = 0; : : : ; M:
The method (5) allows one to nd the items uij recurrently and could be programmed in the form of two nested loops with informational dependance on each other. To nd uij+1; i 2; three values from the previous time layer uij 2; uij 1; uij and two from the current layer uj+1; uij+11 are required.
i 2
For e cient parallel implementation it is necessary to organize the calculation
so that the in uence of the informational dependencies were eliminated. We will
use the approach described in [
We share the task of nding the approximate solution in the grid nodes between processes. Namely, we cut the grid into the K spatial layers. We call the k-th, k = 1; ::; K; spatial layer the set of nodes of the grid as follow f(i; j)k ji = (k
1)N=K + 1; ::; kN=K; j = 1; ::; M g: Let the k-th process calculates a solution in the k-th spatial layer.
The rst process has no informational dependencies on other processes. To nd the solution on the (j + 1)-th time layer the k-th process should receive form the (k 1)-th process two values: u(jk+11)N=K 1 and u(jk+11)N=K ; this constitutes the informational dependencies. Notice, that values u(jk 1)N=K 1 and u(jk 1)N=K j have already been transferred to k-th process and the value u(k 1)N=K+1 have been found.
We call the block Bk;j; k 2; the set of grid nodes f(i; j) ji = (k
1)N=K + 1; ::; kN=K; j = kg: To nd the solution in the block Bk;j; k 2; it is necessary to complete the calculation in the block Bk 1;j: The blocks, where we search the solution in parallel, are on the same separable line. In other words we can iterate through the grid nodes in a direction perpendicular to this line and perform all the operations on the line at each iteration in parallel. Since the task of nding the solution in a single block is complex and data transfer between adjacent processes requires few time, the e ciency of parallelization is high.
Let the K is the number of processes in the pool. In cases of practical importance K m; therefore the situation when the rst process has nished its work and the last process has not started yet is excluded. Let us describe the implementation. In the line 4 the blocking function \receive" is used and non-blocking function \send" in the line 9. for (j from 1 to M)f for (k from 1 to min(j,K))f//To perform in the k-th process if(k>1)
Receive last two dots of Block [k-1,j-k+1]; g g else
Use boundary function; Perform the calculations in the Block [k,j-k+1]; if (k<K)
Send last two dotes of the Block [k,j-k+1] to (k+1)-th process;
The idea of the algorithm is depicted in Fig. 1. Here arrows symbolise the direction in which the approximate solution is sought; little empty rounds represent data, i.e. approximate solution in two nodes, which (k 1)-th process send to k-th one.
To verify the parallel implementation of the algorithm in question we used the test equation + = cos x cos t u(x; t =2) + (x; t); with the following initial and the boundary conditions u(x; t) = sin x cos t; 0 u(0; t) = 0; 0 x t Here the function
is de ned as follow and is identically equal to zero; it was added in the right-hand side of the equation to increase the computational complexity in each iteration.
This initial-boundary problem has the exact solution u(x; t) = sin x cos t: The calculations were performed on a cluster of Ural Branch of RAS \Uran". The following hardware and software con guration were used: CPU INTEL XEON E5450, 2 4 cores, 3GHz; Cache memory 2 6 MB Level 2 cache (5400 Sequence); RAM 16 Gb DDR2; OS Linux 2.6.32; language: C++ with Intel C++ compiler (ICC) v14.0.0; MPI library MVAPICH2 Intel 13.0. We call our computation system architecture as hybrid, since it has distributed memory with shared memory on each computing node.
In the previous paragraph we used the term \process" to refer an abstract unit that performs a calculation; now let us concretize this notion. According to the principles of the number of processes is equal to the number of computational nodes.
Let us present the results of computational experiments for the test equation (see Tab. 1). Recall, that we consider an equidistant partition of [x0; X] and [t0; ]; the grid consists of N and M dotes with respect to space and time. The parameter M was constant through the experiments, M = 4000; parameter N ranges from 500 to 50000 dotes. The number of processes is reported in the rst line and ranges from 1 to 64. Time in seconds is reported in the table cells. Speedup is calculated for multi processes variants (compared to the single process variant).
To characterize the quality of the parallel algorithm the two types of scalability were estimated. Weak scaling is de ned as how the solving time varies with the number of processors for a xed problem size per processor. Weak scaling could be numerically characterized by the function ef fweak(N ) = T (N )=T (N0); where N0 is a xed initial number of processors, in our case N0 = 2, and N is a current number of processors, T (N ) is a time of the solution of the task using N processors. An ideal ef fweak(N ) is identically equal to one and is unreachable according to Amdahl's Law.
Strong scaling is de ned as how the solving time varies with the number of processors for a xed total problem size. Strong scaling could be numerically characterized by the function ef fstrong(N ) = T (N )=(N0 T (N0)): In the case of perfect scalability the computation time is inversely proportional to the number of computing nodes. To estimate the strong scalability the grid of size 50000 5000 was used.
Experimental assessment of function ef fstrong is shown in Fig. 3 (solid line); it is extremely close to perfect scaling (dashed line), which is shown for comparison. 5
For the advection equation with time delay the di erence scheme originally
proposed in [
The parallel version was implemented for the \Uran" supercomputer and could be executed on any number of processors; MPI technology was used. Performed tests proved showed high scaling of the method. We observed the fourteen-fold decrease in time compared with the sequential version. Numerical results are consistent with the theory.
This research is supported by RFBR 14-01-00065 and 13-01-00089. We acknowledge the support by the program 02.A03.21.0006 on 27.08.2013.