# A Computational Realization of a Semi-Lagrangian Method for Solving the Advection Equation

Authors: Alexander Efremov; Eugeniya Karepova; Vladimir Shaydurov; Alexander Vyatkin
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

A parallel implementation of a method of the semi-Lagrangian type for the advection equation on a hybrid architecture computation system is discussed. The difference scheme with variable stencil is constructed on the base of an integral equality between the neighboring time levels. The proposed approach allows one to avoid the Courant-Friedrichs-Lewy restriction on the relation between time step and mesh size. The theoretical results are confirmed by numerical experiments. Performance of a sequential algorithm and several parallel implementations with the OpenMP and CUDA technologies in the C language has been studied.


---

## Body

## 1. Introduction

Many physical phenomena in transport processes are modeled by time-dependent hyperbolic conservation laws [1–4]. The finite volume method (FVM) is a standard conservative method to construct numerical approximations for solving hyperbolic conservation problems. Modern modifications of FVM [5–8] provide well-established conservative methods for solving the governing advection equations. Moreover, some of them were developed to treat high gradients and discontinuities of a solution [7, 8]. In spite of their advances for hyperbolic equations, these methods have the limitation consisting in a time step restriction, mainly for stability sake. On the other hand, during the last three decades the idea of applying the method of characteristics to advective quantities forward in time has rapidly developed and has gained popularity in many areas [9–13]. In contrast to traditional Eulerian schemes, semi-Lagrangian algorithms provide unconditional stability and allow using large time steps. Despite unconditional stability, these methods are explicit and therefore they look well suited for parallelization. Now semi-Lagrangian methods are intensively studied and their efficiency for convection-dominated problems is proved. For a more detailed discussion about the comparison of traditional Eulerian and semi-Lagrangian schemes for hyperbolic conservation laws, see [6, 14, 15].Initially semi-Lagrangian algorithms, as methods of characteristics, were developed with application in climate prediction [16–21]. The simplest schemes use the approximation of a trajectory (or curvilinear characteristic) by a straight line and employ a low-order interpolation to compute a numerical solution. Nowadays, simplicity and efficiency of these schemes make them quite popular in different fields of numerical modeling like fluid dynamics applications [9, 12, 22], shallow water equations [10], fiber dynamics described by the Fokker-Planck equation [11], heat-conduction equation [23], and so forth. Now modern semi-Lagrangian algorithms involve a higher-order approximation of a curvilinear characteristic and employ a higher-order interpolation; see, for example, [22]. Recently, considerable efforts have been made to construct the conservative semi-Lagrangian methods [9, 20, 24–28]. For instance, Scroggs and Semazzi [9] presented a semi-Lagrangian finite volume method that used a rectangular grid for a system of conservation laws and satisfied the discrete conservation relation, but the numerical results demonstrate some violation of full conservation. Early modifications of the semi-Lagrangian approach use a rectangular grid which is fixed throughout the simulation [9, 17, 24, 25]. Semi-Lagrangian schemes allow using spatial grids independently of one another. As a result, adaptive grids are widely used in modern versions of this approach [20, 28]. In spite of the progress in semi-Lagrangian methods, for most of them [9–12, 20, 24, 28] convergence has not been theoretically proved.In this paper, we present a sketch of the theoretical proof and a numerical justification for the difference scheme of the semi-Lagrangian family. We start with the theorem about an exact equality that involves two spatial integrals over domains at neighboring time levels and the third integral over an inflow boundary. To prove convergence theoretically, we use a square grid, the bilinear interpolation, and the Runge-Kutta method for the fourth-order approximation of characteristics. This allows us to prove first-order convergence in a discrete analogue of theL
1-norm. The theoretical convergence estimates are confirmed by numerical results.In the remaining part of the paper, some parallel implementations of this method are studied. We discuss the design subtleties of parallel implementations of the algorithm by the OpenMP-technology for shared memory computational systems and by the CUDA technology for general-purpose GPU programming. In addition, the influence of the HyperThreading technology on the performance of our OpenMP-code is studied. Moreover, the difficulties of the algorithm implementation and the performance for hybrid architecture computation systems are discussed for our CUDA codes.


## 2. Formulation of the Problem

LetD
=
[
0,1
]
×
[
0,1
]. In the closed domain [
0
,
T
]
×
D consider the two-dimensional advection equation
(1)
∂
ρ
∂
t
+
∂
(
u
ρ
)
∂
x
+
∂
(
v
ρ
)
∂
y
=
0
,
where u
(
t
,
x
,
y
) and v
(
t
,
x
,
y
) are known and are sufficiently smooth in [
0
,
T
]
×
D. We suppose for simplicity that ∀
t
∈
[
0
,
T
] the coefficients satisfy the no-slip conditions at the upper and lower sides of D
(2)
u
(
t
,
x
,
y
)
|
y
=
0
=
u
(
t
,
x
,
y
)
|
y
=
1
=
0
,
v
(
t
,
x
,
y
)
|
y
=
0
=
v
(
t
,
x
,
y
)
|
y
=
1
=
0
and the flow conditions at the left and right sides of D
(3)
u
(
t
,
x
,
y
)
|
x
=
0
≥
0
,
u
(
t
,
x
,
y
)
|
x
=
1
≥
0
.
For the unknown function ρ
(
t
,
x
,
y
) the following initial and boundary conditions are specified:
(4)
∀
(
x
,
y
)
∈
D
ρ
(
0
,
x
,
y
)
=
ρ
init
(
x
,
y
)
,
(5)
∀
(
t
,
y
)
∈
[
0
,
T
]
×
[
0,1
]
ρ
(
t
,
0
,
y
)
=
ρ
in
(
t
,
y
)
.


## 3. Numerical Scheme

Subdivide the time segment[
0
,
T
] into K time levels t
k
=
k
τ, k
=
0
,
…
,
K, with the time step τ
=
T
/
K. Let Ω be a closed quadrangle at the time level t
k. For each of its points on the segment t
∈
[
t
k
-
1
,
t
k
] we construct the characteristics defined by the system of ordinary differential equations
(6)
x
′
-
(
t
)
=
u
(
t
,
x
-
(
t
)
,
y
-
(
t
)
)
,
y
′
-
(
t
)
=
v
(
t
,
x
-
(
t
)
,
y
-
(
t
)
)
,
vvvvvvvvv
t
∈
[
t
k
-
1
,
t
k
]
,
with the initial value at level t
=
t
k as a parameter:
(7)
x
-
(
t
k
)
=
x
0
,
y
-
(
t
k
)
=
y
0
,
(
x
0
,
y
0
)
∈
Ω
.
With the help of these characteristics the edges of Ω generate four surfaces S
n, n
=
1
,
…
,
4, with the edges C
n at t
=
t
k
-
1 (Figure 1).Figure 1
Curvilinear quadrangleQ.IfΩ is located near the inflow boundary x
=
0, surfaces S
n can cross the plane x
=
0. In this case we get an additional curvilinear polygon I on the plane x
=
0 (Figure 2).Figure 2
Appearance of the boundary quadrangleI.Generally speaking,I and Q can be triangular, pentagonal, or empty domains. If one of them is empty, then the integral over an empty domain is supposed to be equal to zero. Since there is no fundamental difference, we consider only the most common case with quadrangular domains. For Ω, Q, and I the following statement is valid.Theorem 1.
For a smooth solution of the problem (1)–(5) we have the equality
(8)
∫
Ω
ρ
(
t
k
,
x
,
y
)
d
Ω
=
∫
Q
ρ
(
t
k
-
1
,
x
,
y
)
d
Q
+
∫
I
(
ρ
u
)
(
t
,
0
,
y
)
d
I
.Proof.
Denote the volume bounded byΩ
,
Q
,
I, and surfaces S
n by V and its boundary by Γ. Apply the Gauss-Ostrogradsky theorem to the left-hand side of the equality
(9)
∫
V
(
∂
ρ
∂
t
+
∂
(
u
ρ
)
∂
x
+
∂
(
v
ρ
)
∂
y
)
d
V
=
0
.
Then
(10)
∫
V
(
∂
ρ
∂
t
+
∂
(
u
ρ
)
∂
x
+
∂
(
v
ρ
)
∂
y
)
d
V
=
∫
Γ
(
ρ
,
ρ
u
,
ρ
v
)
·
(
n
t
,
n
x
,
n
y
)
T
d
Γ
=
∫
Γ
ρ
(
1
,
u
,
v
)
·
(
n
t
,
n
x
,
n
y
)
T
d
Γ
=
0
.
Here (
n
t
,
n
x
,
n
y
) is the outer normal to Γ and sign “·” means the scalar product. The normal (
n
t
,
n
x
,
n
y
) equals (
1,0
,
0
) on Ω
,
(
-
1,0
,
0
) on Q, and (
0
,
-
1,0
) on I. For any S
n the normal (
n
t
,
n
x
,
n
y
) is orthogonal to all tangent directions of S
n including the tangent of characteristics (
1
+
u
2
+
v
2
)
-
1
/
2
(
1
,
u
,
v
). Therefore (
1
,
u
,
v
)
·
(
n
t
,
n
x
,
n
y
)
T
=
0. Taking into account this reasoning in (10), we get the equality
(11)
∫
Ω
ρ
(
t
k
,
x
,
y
)
d
Ω
-
∫
Q
ρ
(
t
k
-
1
,
x
,
y
)
d
Q
-
∫
I
(
ρ
u
)
(
t
,
0
,
y
)
d
I
=
0
that is equivalent to the statement of theorem.Now construct the uniform meshD
h with mesh-size h
=
1
/
N, N
≥
2:
(12)
D
h
=
{
(
x
i
,
y
j
)
:
x
i
=
i
h
,
y
j
=
j
h
;
i
,
j
=
0
,
…
,
N
}
.
We will find an approximate solution ρ
h
(
t
,
x
,
y
) at each time level t
=
t
r
∀
r
=
0
,
…
,
K as a grid function with values
(13)
ρ
~
i
,
j
r
=
ρ
h
(
t
r
,
x
i
,
y
j
)
∀
i
,
j
=
0
,
…
,
N
,
unlike the values of the exact solution
(14)
ρ
i
,
j
r
=
ρ
(
t
r
,
x
i
,
y
j
)
.To construct the difference scheme with a variable stencil, we suppose that the functionρ
h at time level t
k
-
1 is already known and we need to find it at level t
k. To compute ρ
~
i
,
j
k for some i
,
j
=
1
,
2
,
…
,
N
-
1, we take the square Ω
i
,
j with four vertices (
x
i
±
h
/
2
,
y
j
±
h
/
2
) and apply Theorem 1. To determine ρ
~
i
,
j
k on the boundary of D we use the rectangles Ω
i
,
j which are adjoined to this boundary inside D (Figure 3).Figure 3
Boundary rectangles.Note thatρ
~
0
,
j
k
=
ρ
0
,
j
k are known from the boundary condition (5).Without loss of generality we describe the construction of the difference equations for inner nodes withi
,
j
=
1
,
2
,
…
,
N
-
1 only. Thus, due to Theorem 1 we get
(15)
∫
Ω
i
,
j
ρ
(
t
k
,
x
,
y
)
d
Ω
=
∫
Q
i
,
j
k
-
1
ρ
(
t
k
-
1
,
x
,
y
)
d
Q
+
∫
I
i
,
j
k
-
1
(
ρ
u
)
(
t
,
0
,
y
)
d
I
.
Here the curvilinear polygons Q
i
,
j
k
-
1 and I
i
,
j
k
-
1 are formed by the characteristics (6) that issue out of the edges of square Ω
i
,
j. To compute the second integrals in (15) numerically, first we replace the exact function ρ
(
t
k
-
1
,
x
,
y
) by its bilinear interpolant
(16)
ρ
h
I
(
t
k
-
1
,
x
,
y
)
=
∑
q
=
0
N
∑
p
=
0
N
ρ
p
,
q
k
-
1
φ
p
,
q
(
x
,
y
)
with the help of the basis functions
(17)
φ
p
,
q
(
x
)
=
{
(
1
-
|
x
p
-
x
|
h
)
(
1
-
|
y
q
-
y
|
h
)
∀
(
x
,
y
)
∈
[
x
p
-
1
,
x
p
+
1
]
×
[
y
q
-
1
,
y
q
+
1
]
,
0
otherwise
,
∀
p
,
q
=
0
,
…
,
N
.
vv
To compute the integral over the domain I
i
,
j
k
-
1, we also use the bilinear interpolant
(18)
(
ρ
u
)
τ
I
(
t
,
y
)
=
∑
q
=
0
N
∑
r
=
0
K
ρ
in
(
t
r
,
y
q
)
u
(
t
r
,
0
,
y
q
)
ψ
r
,
q
(
t
,
y
)
with the basis functions
(19)
ψ
r
,
q
(
x
)
=
{
(
1
-
|
t
r
-
t
|
τ
)
(
1
-
|
y
q
-
y
|
h
)
∀
(
t
,
y
)
∈
[
t
r
-
1
,
t
r
+
1
]
×
[
y
q
-
1
,
y
q
+
1
]
,
0
otherwise
,
∀
r
=
0
,
…
,
K
,
q
=
0
,
…
,
N
.
vv
The left-hand side of (15) is approximated by the midpoint quadrature rule with second-order accuracy:
(20)
∫
Ω
i
,
j
ρ
(
t
k
,
x
,
y
)
d
Ω
≈
mes
(
Ω
i
,
j
)
ρ
i
,
j
k
∀
i
,
j
=
0,1
,
…
,
N
.
So, instead of the exact equality (15) we get the approximate one:
(21)
mes
(
Ω
i
,
j
)
ρ
i
,
j
k
≈
∫
Q
i
,
j
k
-
1
ρ
h
I
(
t
k
-
1
,
x
,
y
)
d
Q
+
∫
I
i
,
j
k
-
1
(
ρ
u
)
τ
I
(
t
,
0
,
y
)
d
I
.
To simplify the numerical computation of the right-hand side in (21), we approximate the domains Q
i
,
j
k
-
1 and I
i
,
j
k
-
1 by simpler ones. Since in the more general case both domains are curvilinear quadrangles, we demonstrate the approximation only for quadrangular Q
i
,
j
k
-
1. Introduce four additional points (
x
i
±
h
/
2
,
y
j
) and (
x
i
,
y
j
±
h
/
2
) on the square Ω
i
,
j at time level t
k and denote each of the eight nodes by A
n
=
(
x
^
n
,
y
^
n
), n
=
1
,
…
,
8. Out of each A
n we pass the corresponding characteristic to the time level t
k
-
1 which gives the point B
n
=
(
x
-
n
,
y
-
n
) (Figure 4).Figure 4
Approximation of a curvilinear quadrangle.To compute the coordinates of the pointB
n numerically, we solve the system of ordinary differential equations (6) with the initial condition
(22)
x
~
(
t
k
)
=
x
n
,
y
~
(
t
k
)
=
y
n
,
by the fourth-order Runge-Kutta method [29]. Thus, we find the approximation B
n
h
=
(
x
~
n
(
t
k
-
1
)
,
y
~
n
(
t
k
-
1
)
) of the point B
n. The nodes B
n
h, n
=
1
,
…
,
8, define the polygon P
i
,
j
k
-
1 which is considered as a quadrangle with four parabolic edges (Figures 4-5). The constructed domain P
i
,
j
k
-
1 approximates Q
i
,
j
k
-
1. In the same way we construct the polygon L
i
,
j
k
-
1 which approximates I
i
,
j
k
-
1. For the above approximation the following statement is valid [30].Figure 5
Approximation of nodes and edges.Lemma 2.
Let the coordinates of the nodesB
n, n
=
1
,
…
,
8, be computed within the fourth-order accuracy x
-
n
-
x
~
n
(
t
k
)
=
O
(
τ
4
),  y
-
n
-
y
~
n
(
t
)
=
O
(
τ
4
). Assume that the ratio between τ and h is fixed: τ
=
c
~
h. Then for all i
,
j
=
0
,
…
,
N
(23)
mes
(
Q
i
,
j
k
-
1
∖
P
i
,
j
k
-
1
)
+
mes
(
P
i
,
j
k
-
1
∖
Q
i
,
j
k
-
1
)
=
O
(
h
4
)
,
mes
(
I
i
,
j
k
-
1
∖
L
i
,
j
k
-
1
)
+
mes
(
L
i
,
j
k
-
1
∖
I
i
,
j
k
-
1
)
=
O
(
h
4
)
,
where the notation mes
(
Ω
) means the measure of the domain Ω.Thus, the replacement ofQ
i
,
j
k
-
1 by P
i
,
j
k
-
1 and I
i
,
j
k
-
1 by L
i
,
j
k
-
1, i
=
1
,
…
,
N, j
=
0
,
…
,
N, reduces the approximate equality (21) to another one:
(24)
mes
(
Ω
i
,
j
)
ρ
i
,
j
k
≈
∫
P
i
,
j
k
-
1
ρ
h
I
(
t
k
-
1
,
x
,
y
)
d
P
+
∫
L
i
,
j
k
-
1
(
ρ
u
)
τ
I
(
t
,
0
,
y
)
d
L
.
Divide it by mes
(
Ω
i
,
j
) and replace ρ
h
I by the interpolant ρ
~
h
I of the known grid function ρ
h at the level t
k
-
1:
(25)
ρ
~
h
I
(
t
k
-
1
,
x
,
y
)
=
∑
q
=
0
N
∑
p
=
0
N
ρ
~
p
,
q
k
-
1
φ
p
,
q
(
x
,
y
)
.
As a result, we get the equation for finding ρ
~
i
,
j
k as an approximation of ρ
i
,
j
k:
(26)
ρ
~
i
,
j
k
=
1
mes
(
Ω
i
,
j
)
∫
P
i
,
j
k
-
1
ρ
~
h
I
(
t
k
-
1
,
x
,
y
)
d
P
+
1
mes
(
Ω
i
,
j
)
∫
L
i
,
j
k
-
1
(
ρ
u
)
τ
I
(
t
,
0
,
y
)
d
L
vvvvvv
∀
i
=
1
,
…
,
N
,
j
=
0,1
,
…
,
N
,
(27)
ρ
~
0
,
j
k
=
ρ
in
(
t
k
,
0
,
y
j
)
∀
j
=
0
,
…
,
N
.
To compute the integrals numerically, we decompose the domain P
i
,
j
k
-
1 (or L
i
,
j
k
-
1) into several triangles which lie only in one cell [
x
p
-
1
,
x
p
]
×
[
y
q
,
y
q
+
1
], p
,
q
=
0,1
,
…
,
N
-
1, and have only one parabolic edge and two straight-line edges parallel to the coordinate axes. Then we replace the integral over the domain P
i
,
j
k
-
1 (or L
i
,
j
k
-
1) by a sum of integrals. Thus we compute integrals directly without any quadrature rule.To evaluate the order of convergence, we use the discrete analogue of theL
1
(
D
)-norm:
(28)
∥
ρ
h
∥
L
1
h
=
∑
i
,
j
=
0
N
|
ρ
i
,
j
h
|
mes
(
Ω
i
,
j
)
.
For the numerical solution computed by (26) the following theorem is valid.Theorem 3.
Let the solutionρ
(
t
,
x
,
y
) of the problem (1)–(5) be sufficiently smooth and let the discrete solution ρ
h be computed by (26). Assume that τ
=
c
~
h. Then we have the following estimate: ∀
k
=
0
,
1
,
…
,
K
(29)
∥
ρ
(
t
k
,
·
)
-
ρ
h
(
t
k
,
·
)
∥
L
1
h
≤
k
(
c
1
h
2
+
c
2
c
~
2
h
3
)
,
with the constants c
1 and c
2 independent of k, h, τ, and c
~.Proof.
Use the induction onk. For k
=
0 inequality (29) is valid because of the exact initial condition (4). Suppose the estimate (29) is valid for some k
-
1
≥
0 and prove it for k.
From Theorem1  ∀
i
=
1
,
…
,
N, j
=
0
,
…
,
N we get
(30)
∫
Ω
i
,
j
ρ
(
t
k
,
x
,
y
)
d
Ω
=
∫
Q
i
,
j
k
-
1
ρ
(
t
k
-
1
,
x
,
y
)
d
Q
+
∫
I
i
,
j
k
-
1
ρ
(
t
,
0
,
y
)
u
(
t
,
0
,
y
)
d
I
.
For the first integral we have the equality
(31)
ρ
(
t
k
,
x
i
,
y
j
)
mes
(
Ω
i
,
j
)
=
∫
Ω
i
,
j
ρ
(
t
k
,
x
,
y
)
d
Ω
+
δ
i
,
j
k
mes
(
Ω
i
,
j
)
,
where
(32)
|
δ
i
,
j
k
|
≤
{
c
~
1
h
2
∀
i
=
1
,
…
,
N
-
1
,
c
~
2
h
when
i
=
N
.
In the second and third integrals we replace the polygonsQ
i
,
j
k
-
1 and I
i
,
j
k
-
1 by P
i
,
j
k
-
1 and L
i
,
j
k
-
1, respectively, according to the foregoing approximation. Then we replace ρ
(
t
k
-
1
,
x
,
y
) and ρ
(
t
,
0
,
y
)
u
(
t
,
0
,
y
) by their piecewise bilinear interpolants. Due to Lemma 2 and boundedness of the functions ρ
(
t
,
x
,
y
) and u
(
t
,
x
,
y
), formulae (30)-(31) are represented in the following way:
(33)
ρ
i
,
j
k
mes
(
Ω
i
,
j
)
=
∫
P
i
,
j
k
-
1
ρ
h
I
(
t
k
-
1
,
x
,
y
)
d
P
+
∫
L
i
,
j
k
-
1
(
ρ
u
)
τ
I
(
t
,
0
,
y
)
d
L
+
δ
i
,
j
k
mes
(
Ω
i
,
j
)
+
γ
i
,
j
k
-
1
+
η
i
,
j
k
-
1
mes
(
P
i
,
j
k
-
1
)
+
θ
i
,
j
k
-
1
mes
(
L
i
,
j
k
-
1
)
,
where |
γ
i
,
j
k
-
1
|
≤
c
~
3
h
4, |
η
i
,
j
k
-
1
|
≤
c
~
4
h
2, |
θ
i
,
j
k
-
1
|
≤
c
~
5
τ
h. Now multiply (26) by mes
(
Ω
i
,
j
) and subtract it from (33). Then we have
(34)
(
ρ
i
,
j
k
-
ρ
~
i
,
j
k
)
mes
(
Ω
i
,
j
)
=
∫
P
i
,
j
k
-
1
∑
p
,
q
=
0
N
(
ρ
p
,
q
k
-
1
-
ρ
~
p
,
q
k
-
1
)
ψ
p
,
q
(
x
,
y
)
d
P
+
δ
i
,
j
k
mes
(
Ω
i
,
j
)
+
γ
i
,
j
k
-
1
+
η
i
,
j
k
-
1
mes
(
P
i
,
j
k
-
1
)
+
θ
i
,
j
k
-
1
mes
(
L
i
,
j
k
-
1
)
.
Now, let us sum up the modulus of the last equation for all i
=
1
,
…
,
N, j
=
0
,
…
,
N and use the decomposition
(35)
ρ
~
i
,
j
k
-
1
=
ρ
i
,
j
k
-
1
+
ξ
i
,
j
k
-
1
at level t
k
-
1 with a grid function ξ
k
-
1
(
x
,
y
) satisfying the estimate
(36)
∥
ξ
k
-
1
∥
L
1
h
≤
(
k
-
1
)
(
c
1
h
2
+
c
2
c
~
2
h
3
)
,
due to the induction hypothesis. Then we get
(37)
∥
ξ
k
∥
L
1
h
≤
∑
i
,
j
=
0
N
(
∫
P
i
,
j
k
-
1
∑
p
,
q
=
0
N
|
ξ
p
,
q
k
-
1
|
ψ
p
,
q
(
x
,
y
)
d
P
)
+
∑
i
,
j
=
0
N
|
δ
i
,
j
k
|
mes
(
Ω
i
,
j
)
+
|
γ
i
,
j
k
-
1
|
+
∑
i
,
j
=
0
N
(
|
η
i
,
j
k
-
1
|
mes
(
P
i
,
j
k
-
1
)
+
|
θ
i
,
j
k
-
1
|
mes
(
L
i
,
j
k
-
1
)
)
.
Since
(38)
∑
i
,
j
=
0
N
|
δ
i
,
j
k
|
mes
(
Ω
i
,
j
)
≤
(
c
~
1
+
2
c
~
2
)
h
2
,
∑
i
,
j
=
0
N
mes
(
P
i
,
j
k
-
1
)
≤
1
,
∑
i
,
j
=
0
N
mes
(
L
i
,
j
k
-
1
)
=
τ
,
we have
(39)
∥
ξ
k
∥
L
1
h
≤
∑
i
,
j
=
0
N
(
∫
P
i
,
j
k
-
1
∑
p
,
q
=
0
N
|
ξ
p
,
q
k
-
1
|
ψ
p
,
q
(
x
,
y
)
d
P
)
+
(
c
~
1
+
2
c
~
2
+
c
~
3
+
c
~
4
)
h
2
+
c
~
5
h
τ
2
.
Finally consider the transformations
(40)
∑
i
,
j
=
0
N
(
∫
P
i
,
j
k
-
1
∑
p
,
q
=
0
N
|
ξ
p
,
q
k
-
1
|
ψ
p
,
q
(
x
,
y
)
d
P
)
=
∑
p
,
q
=
0
N
(
∑
i
,
j
=
0
N
∫
P
i
,
j
k
-
1
|
ξ
p
,
q
k
-
1
|
ψ
p
,
q
(
x
,
y
)
d
P
)
≤
∑
p
,
q
=
0
N
(
|
ξ
p
,
q
k
-
1
|
∫
D
ψ
p
,
q
(
x
,
y
)
d
D
)
.
It leads to the inequality
(41)
∥
ξ
k
∥
L
1
h
≤
∑
p
,
q
=
0
N
|
ξ
p
,
q
k
-
1
|
mes
(
Ω
p
,
q
)
+
(
c
~
1
+
2
c
~
2
+
c
~
3
+
c
~
4
)
h
2
+
c
~
5
h
τ
2
.
Denote c
1
=
(
c
~
1
+
2
c
~
2
+
c
~
3
+
c
~
4
) and c
2
=
c
~
5. Thus due to the relation τ
=
c
~
h we get inequality (29).Corollary 4.
Let the conditions of Theorem3 be valid. Then for t
k
=
T one has the estimate
(42)
∥
ρ
(
T
,
·
)
-
ρ
h
(
T
,
·
)
∥
L
1
h
≤
T
(
c
1
h
c
~
+
c
2
c
~
h
2
)
.Remark 5.
Let the functionsρ
init
(
x
,
y
) and ρ
in
(
t
,
y
) be greater than zero in the initial and boundary conditions (4)-(5). Then the interpolants ρ
~
h
I
(
t
0
,
x
,
y
)  and  (
ρ
u
)
τ
I
(
t
,
0
,
y
) due to (3) are nonnegative. The integration of them results in nonnegative values in (26). Thus, by induction we can prove the inequality
(43)
ρ
h
(
t
r
,
x
i
,
y
j
)
≥
0
∀
r
=
1
,
…
,
K
,
i
,
j
=
0
,
…
,
N
.Remark 6.
The strategy of the domain approximation with 8 nodes is not optimal, of course. In fact it is enough to use 4 nodes for rectangles. But a theoretical justification becomes much more complicated. We did not demonstrate it here since the main purpose of the paper is to verify parallel properties of the proposed algorithm. Of course, such an optimization reduces the number of arithmetical operations but has no influence on the parallel properties of the algorithm. The same applies to the difference schemes of higher order.


## 4. The Numerical Algorithm and Its Parallel Implementations

The constructed algorithm is implemented in the following way.Algorithm 7 (sequential).
(1)
Setρ
h
(
0
,
x
i
,
y
j
)
=
ρ
init
(
x
i
,
y
j
)
,
i
,
j
=
0
,
…
,
N, as the initial data (4).
(2)
Time loop: for each time step k
=
1
,
…
,
K do:
Space loop: for each cell Ω
i
,
j
,
i
=
1
,
…
,
N
,
j
=
0
,
…
,
N do:
(2.1)
For each nodeA
n
=
(
x
^
n
,
y
^
n
)
,
n
=
1
,
…
,
4, solve the system (18)–(20) and determine the corresponding vertex coordinates B
n
h
=
(
x
~
n
(
t
k
-
1
)
,
y
~
n
(
t
k
-
1
)
) of a polygon P
i
,
j
k
-
1.
(2.2)
If a certain characteristicA
n
B
n
h intersects the plane x
=
0 then the coordinates of this cross-point are determined.
(2.3)
Compute(44)
J
i
,
j
k
-
1
=
∫
P
i
,
j
k
-
1
ρ
h
I
(
t
k
-
1
,
x
,
y
)
d
P
+
∫
L
i
,
j
k
-
1
(
ρ
u
)
I
(
t
,
0
,
y
)
d
L
,
according to (26) where the integrals are calculated over each nonempty intersection P
i
,
j
k
-
1
∩
{
[
x
p
,
x
p
+
1
]
×
[
y
q
,
y
q
+
1
]
} and L
i
,
j
k
-
1
∩
{
[
t
k
-
1
,
t
k
]
×
[
y
q
,
y
q
+
1
]
} separately.
(2.4)
Computeρ
h
(
t
k
,
x
i
,
y
j
)
=
J
i
,
j
k
-
1
/
mes
(
Ω
i
,
j
).
The end of the space loop.
(2.5)
Putρ
h
(
t
k
,
0
,
y
j
)
=
ρ
in
(
t
k
,
y
j
)  for  all  j
=
0
,
…
,
N.
(2.6)
If necessary, calculate the norms of a solution, an error, and other statistic data with respect to an actual time step.
The end of the time loop.Note that items (2.1)–(2.3) are compute-intensive, especially the procedure of determining the mutual arrangement ofP
i
,
j
k
-
1 and the cells {
[
x
p
,
x
p
+
1
]
×
[
y
q
,
y
q
+
1
]
}
p
,
q
=
0
N
-
1 at the previous time level in item (2.3).The algorithm is explicit with respect to time, since to calculateρ
h
(
t
k
,
x
,
y
) at each time level t
k the data are used only from the previous time level t
k
-
1.Another advantage of Algorithm7 is data independence in the general space loop; that is, the items (2.1)–(2.4) are carried out for any pair (
i
,
j
)
,
i
=
1
,
…
,
N
,
j
=
0
,
…
,
N, independently. In this connection the data parallelism is used.In the shared memory case for the OpenMP-technology it is sufficient to parallelize the general space loop at each time level using an OpenMP directive like the following one:#pragma omp parallel for collapse (2)…In order for paralleling to be correct, the data-sharing attributes of all variables to intermediate outcomes of items (2.1)–(2.4) have to be private for each thread.Another justified approach to paralleling the algorithm is the usage of the NVIDIA CUDA technology for GPU. The main aspects of the parallel implementation related to various features of general-purpose GPU programming are briefly discussed below.All functions used in the numerical calculations on a CPU must be recompiled for a GPU. If we use NVIDIA CUDA these functions must be declared with the special qualifier  __host__ __device__ which indicates that the NVCC compiler creates two versions of its executable code for a CPU (host) and for a GPU (device) separately. GPU will call a device version of a function while CPU will call its host version.The principles of efficient CUDA programming are as follows: (1) the maximal use of inherent parallelism of the problem and (2) the optimization of memory access.The first version of our parallel CUDA-algorithm is based on the inherent parallelism of our numerical explicit approach. Every thread treats only one cellΩ
i
,
j
∈
D
h; hence, the space loop body (items (2.1)–(2.4) of Algorithm 7) is the general computation kernel.While programming for GPU, the correct defining of the kernel configuration is important. The kernel configuration includes two parameters, namely, the number of blocks (blockCount) in a grid and the number of threads (blockSize) per block. There is a limit in 1024 threads per block for our GPU NVIDIA hardware. In the first CUDA-algorithm no threads amount optimization is used.Consequently, the simplest parallel version of Algorithm7 for the CPU/GPU hybrid architecture is the following.Algorithm 8 (CUDA parallel, version 1).
(1)
Calculate the kernel configuration (blockSize,blockCount) using data about a computational domain.
(2)
Allocate host (CPU) and device (GPU) memory; copy initial data from host to device.
(3)
Time loop: for each time step k
=
1
,
…
,
K do:
(3.1)
Call the first CUDA kernel (basic):
(3.1.1)
For each cellΩ
i
,
j,  i
=
1
,
…
,
N,  j
=
0
,
…
,
N, execute items (2.1)–(2.4) of Algorithm 7 in parallel.
(3.2)
Synch point: wait for calculations to be completed.
(3.3)
Call the second CUDA kernel (assistive):
(3.3.1)
Copy data from an actual time level array to the previous one in parallel.
(3.4)
Synch point: wait for copying to be completed.
(3.5)
If necessary, copy results from device to host.
(3.6)
If necessary, calculate the norms of a solution, an error, and other statistic data with respect to an actual time step.
The end of the time loop.
(4)
Copy the results from device to host.In order to decrease the execution time of Algorithm8, items (3.5)-(3.6) must be performed as rare as possible, for instance, only at time levels where accuracy control, data for drawing, and so forth are needed.Algorithm8 has two general disadvantages: (1) small speedup in comparison to the sequential version (Figure 8) and (2) impossibility of execution with a fine computational mesh (Table 2). What is the matter with these problems?First, the general loop has a lot of selection statements.The main selection is between two different ways of catching in item (3.1.1) of Algorithm8 (to be more exact, in the item (2.3) of sequential Algorithm 7). The cells whose trajectories intersect the boundary and the internal cells are processed in different ways. We can use two different kernels which execute in parallel.We use data parallelism only in Algorithm8. However, there is yet another class of parallelism to be exploited on NVIDIA GPU. This parallelism is similar to the task parallelism that is found in multithreaded CPU applications. NVIDIA CUDA task parallelism is based on CUDA streams. A CUDA stream represents a queue of GPU operations such as kernel launches, memory copies, and event starts and stops. The order in which operations are added to the stream specifies the order in which they will be executed. Each stream may be considered as a certain task on the GPU, and there are opportunities for these tasks to execute in parallel [31]. Thus we apply CUDA streams to our two kernels to improve parallelism and total GPU utilization.There are many selections in item (2.3) for determining mutual arrangement ofP
i
,
j
k
-
1 and cells of the mesh of the previous time level. Unfortunately, we cannot avoid these selections.Secondly, the CUDA kernel in (3.3) of Algorithm8 idles in the context of computation. We apply loop unrolling in order to eliminate this kernel.Thirdly, the basic CUDA kernel of Algorithm8 is not optimal for memory access. To optimize concurrent read access global memory of simultaneously running threads, constant memory is preferable to use. Applying this approach in our program we allocate all invariable values in constant memory GPU.Moreover, for the optimization of parameters of kernels launch we use the CUDA Occupancy Calculator. It calculates optimal streaming multiprocessor (SM) utilization taking into account GPU compute capability, CUDA device properties, the number of blocks in a grid, the number of threads per block, the size of the shared-memory space, and the number of registers per thread.Consequently, the second CUDA-version of Algorithm7 for the CPU/GPU hybrid architecture is the following.Algorithm 9 (CUDA parallel, version 2).
(1)
Calculate kernel configuration (blockSize,blockCount) using data about a computational domain.
(2)
Allocate host (CPU) and device (GPU) memory, copy initial data from host to device, and copy constant data from host to constant memory of device.
(3)
Time loop: for each time step k
=
1,3
,
5
,
…
,
K
-
1 do:
(3.1)
Call the first CUDA kernel to boundary cells access by the first CUDA stream:
(3.1.1)
Execute items (2.1)–(2.4) of Algorithm7 in parallel for each cell Ω
i
,
j whose characteristics intersect the boundary.
(3.2)
Call the second CUDA kernel to inner cells access by the second CUDA stream:
(3.2.1)
Execute items (2.1), (2.3)-(2.4) of Algorithm7 in parallel for each internal cell Ω
i
,
j.
(3.3)
Synch point: wait for calculations of both kernels to be completed.
(3.4)
Call the first CUDA kernel to boundary cells access by the first CUDA stream:
(3.4.1)
Execute items (2.1)–(2.4) of Algorithm7 in parallel for each cell Ω
i
,
j whose characteristics intersect the boundary.
(3.5)
Call the second CUDA kernel to inner cells access by the second CUDA stream:
(3.5.1)
Execute items (2.1), (2.3)-(2.4) of Algorithm7 in parallel for each internal cell Ω
i
,
j.
(3.6)
Synch point: wait for calculations of both kernels to be completed.
(3.7)
If necessary, copy results from device to host.
(3.8)
If necessary, calculate the norms of a solution, an error, and other statistic data with respect to an actual time step.
The end of the time loop.
(4)
IfK is odd then repeat items (3.1)-(3.2).
(5)
Copy results from device to host.


## 5. Numerical Experiments

Specify the velocities(45)
u
(
t
,
x
,
y
)
=
100
y
(
1
-
y
)
[
π
2
-
arctg
(
x
)
]
,
v
(
t
,
x
,
y
)
=
arctg
(
x
(
1
-
x
)
y
(
1
-
y
)
(
1
+
t
)
10
)
and take the initial and boundary conditions in the following form:
(46)
∀
(
x
,
y
)
∈
D
ρ
(
0
,
x
,
y
)
=
ρ
init
(
x
,
y
)
=
1
,
∀
(
t
,
y
)
∈
[
0
,
T
]
×
[
0,1
]
ρ
(
t
,
0
,
y
)
=
ρ
in
(
t
,
y
)
=
1
.Numerical experiments were performed with the ICM SB RAS FLAGMAN computation system of the following configuration.Hardware. CPU: INTEL XEON, 2
×
6 cores, HyperThreading; 40 Gb DDR3; GPU: NVIDIA TESLA C2050 (CC 2.0).Software. OS: UBUNTU 11.04; C/C++: GCC 4.5.2, INTEL C++ Compiler 13.1.0; CUDA C/C++: NVCC 5.0; NVIDIA CUDA 5.0; BOOST 1.53; NVIDIA CUDA-GDB 5.0.One of the purposes for the numerical experiments was to check the convergence order inτ and h. Therefore the computations were performed on the sequence of N
×
N regular square grids, N
=
10
·
2
n, n
=
0
,
…
,
6. The number of time steps is defined by τ
=
h
/
5.Assume that{
ρ
h
n
}
n
=
0
6 is the set of solutions found on the sequence of square grids. The expression lo
g
2
(
∥
ρ
-
ρ
h
n
∥
L
1
h
/
∥
ρ
-
ρ
h
n
+
1
∥
L
1
h
) as a function in n can be considered as the order of convergence (Figure 6). The corresponding exact values of ρ
i
,
j
K were computed by the characteristics method directly. In Figure 6 we can see the first order of convergence in h.Figure 6
The order of convergence.In our sequential and OpenMP computational experiments we compare the GCC and Intel C++ compilers. The execution time for the better Intel compiling code is on average by 15% less than the better GCC one. All presented numerical results were compiled with −O2 optimization level. Let us remark that we try to use other compiler options (−O3, −parallel, −AVX), but the performance increases only slightly or sometimes even decreases. For the CUDA-version we used the NVCC compiler.The results of computation speedup of the OpenMP-version are presented in Table1 and in Figure 8. The first line of the table shows speedup (or rather slowing down) of one thread that executed the OpenMP-code as compared with the sequential code. Generally, the compiler optimization under the OpenMP library and overhead related to the synchronization of OpenMP-threads can be estimated for these data. As we can see in our case, overhead is small.Table 1
Speedup of OpenMP-code (HyperThreading is switched Off/On).
Number of threads
Number of points in one space dimension (number of the time steps)
80
*
80 (400)
160
*
160 (800)
320
*
320 (1600)
640
*
640 (3200)
1280
*
1280 (6400)
1
0.94/1.00
0.99/1.00
0.99/1.00
0.99/1.00
0.99/1.00
4
3.86/3.89
3.90/3.93
3.92/3.93
3.93/3.95
3.93/3.96
8
7.40/7.35
7.48/7.37
7.53/7.57
7.55/7.59
7.56/7.60
12
10.86/10.86
11.08/11.04
11.17/11.12
11.21/11.14
11.20/10.86
16
4.15/8.76
4.77/8.94
6.32/8.96
7.54/8.98
7.60/8.88
20
4.69/10.63
5.82/10.85
6.55/10.94
7.26/10.96
7.89/10.73
24
5.71/10.98
6.17/12.80
7.30/12.55
7.97/12.83
8.45/12.85Table 2
Execution time of all versions of program.
Version
Number of mesh points in one space dimension (number of the time steps)
80
*
80 (400)
160
*
160 (800)
320
*
320 (1600)
640
*
640 (3200)
1280
*
1280 (6400)
Sequential, −O0
20.16
159.40
1268.46
10107.80
*
Sequential, −O2
9.99
78.97
626.72
4980.61
39598.90
Sequential, −O3
9.87
78.08
619.80
4936.25
39202.91
OpenMP(12), −O0, HT Off
1.98
12.52
103.45
819.06
6519.87
OpenMP(12), −O2, HT Off
0.92
7.13
56.10
444.15
3535.53
OpenMP(24), −O2, HT On
0.91
6.17
49.93
388.20
3080.91
CUDA version 1
2.55
13.06
74.02
**
**
CUDA version 2
3.20
8.27
42.60
308.82
**One of the purposes of the studies is to assess influence of the HyperThreading (HT) technology on the parallelism. As long as only 12 physical cores are available on the CPU, we have access to 24 logical cores when HT is enabled and to 12 logical cores when HT is disabled, respectively. Experiments show that when HT is enabled the execution time with 24 threads is about 14% less than that with 12 threads when HT is disabled (Figure7). As our code is compute-intensive, probably, the advantage of HT may be related to optimization of memory access.Figure 7
Execution time of OpenMP-code (the main vertical axis) and speedup in comparison with the sequential code (additional vertical axis). The comparison of results when HyperThreading is switched On or Off.Figure 8
Speedup of parallel versions.The execution time of the sequential, OpenMP, and two versions of CUDA codes is given in Table2. The ≪
*
≫ symbol signifies that the result has not been reached in reasonable time; for the CUDA-versions the ≪
*
*
≫ symbol means that the kernel launch failed due to registers bottleneck. For OpenMP the results on 12 threads when HT is disabled and on 24 threads when HT is enabled are demonstrated.Comparative information on a possibility of code optimization by the GCC compiler is also presented in Table2. The execution time of the code compiled without optimization and with −O2 and −O3 optimization levels, respectively, was measured. It is clear from Table 2 that compiler optimization considerably (more than two times) reduces the execution time of the sequential code. As the compiler does not optimize the CUDA-code, the execution time of the CUDA-versions does not depend on compiling options.Table3 and Figure 8 present the data on computation speedup of the best OpenMP-version and two CUDA-versions in comparison with the sequential program compiled with −O2. Speedup of the second CUDA-version in comparison with the OpenMP-version is given as well. Numerical experiments show that on fine grids in comparison with the sequential program the best Open-MP program with 24 threads gives more than 12 times speedup and the second CUDA-version shows about 16 times speedup.Table 3
Speedup of parallel versions.
Number of mesh points in one space dimension  (number of the time steps)
80
*
80 (400)
160
*
160 (800)
320
*
320 (1600)
640
*
640 (3200)
OpenMP(24), -O2, HT ON/sequential
10.98
12.80
12.55
12.83
CUDA version 1/sequential
3.92
6.05
8.47
**
CUDA version 2/sequential
3.12
9.55
14.71
16.13
CUDA version 2/OpenMP(12), -O2, HT Off
0.29
0.86
1.32
1.44
CUDA version 2/OpenMP(24), -O2, HT On
0.28
0.75
1.17
1.26


## 6. Discussion

Nowadays there are a lot of algorithms of the family of semi-Lagrangian methods. As we mentioned above, the presented method is based on square grid only and uses accessory algorithms that makes computation more resource-intensive. Nevertheless, such a complication allows us to prove first-order convergence. Furthermore, the theorem which allows one to take into account a volume of substance passed through a boundary is presented. This enables us to prove the balance equation. Numerical experiments completely confirm the theoretical convergence results.As for parallel implementation of our approach, we can note the following.Though the algorithm is explicit we are not satisfied with the results of the CUDA-versions.First of all, there is a problem with feasibility of CUDA-code on fine grids (more than 640 × 640 nodes). Profile-feedback analysis shows at least two causes: (1) assumed computational domain decomposition and (2) the problem of hardware constraint on the amount of registers on streaming multiprocessors.Concerning the first item, it should be noted that in our approach a 2D computational domain is mapped to a 1D array of cells, and every thread treats one cell (e.g., see Pseudocode1).Pseudocode 1
__global__  void  unit_foo  (float*  inA,  float  out,  int  size)  {
int  i=blockIdx.x*blockDim.x+threadIdx.x;
if  (i<size)
out[i]  =  do_smb_with(inA[i]);
}For the parameters of kernels launch to be readily adaptable, we can apply “thread reuse”; namely, every thread treats some uncoupled cells (see Pseudocode2).Pseudocode 2
__global__  void  vec_foo  (float*  inA,  float  out,  int  size)  {
for  (int  i=blockIdx.x*blockDim.x+threadIdx.x;
i<size;  i+=blockDim.x*gridDim.x)
out[i]  =  do_smb_with(inA[i]);
}Besides, we can use several video adapters. In this case we should resolve the problem of the distribution of computational cells between devices under the following restriction: we do not know in advance how many cells of the previous time level are required to calculate an actual value (e.g., varying-width shadow line).The problem with registers in our case is related to deep nesting level of functions calls in the item (2.3) for determining the mutual arrangement ofP
i
,
j
k
-
1 and cells of a mesh of the previous time level. Indeed, this is a bottleneck of our sequential algorithm. To resolve this issue, we should modify the initial sequential algorithm, unfortunately.In addition, the item (2.3) has many flow control instructions (“if” statements, mainly), but we cannot say that this branching creates some significant performance penalty in terms of SIMT architecture. Indeed, in a CUDA kernel, any flow control instruction (if, switch, do, for, while) can significantly affect the instruction throughput by causing threads of the same warp to diverge [32]. If this happens, the different execution paths must be serialized, since all of the threads of a warp share a program counter; this increases the total number of instructions executed for this warp. When all the different execution paths have completed, the threads converge back to the same execution path.Fortunately, this rule works in the cases where the control flow depends on the thread ID only. However, in our case we have many branches which do not depend on thread ID. Inside our computational kernel, the main part of branches looks like Pseudocode3.Pseudocode 3
if  ((indCurSqOx[1]  >=  0)  &&  (indCurSqOy[1]  >=  0))  {
do_smth();
}else  {
do_smth_else();
}As we can see, our “IF” statement does not depend on a thread ID. Thus, this type of branch does not affect the performance of the CUDA kernel.


## 7. Conclusion

We present an unconditionally stable semi-Lagrangian scheme of the first-order accuracy. The numerical experiments confirm the theoretical studies. Performance of sequential and several parallel algorithms implemented with the OpenMP and CUDA technologies in the C language was studied. The optimization potential of the CUDA-version remains unexhausted yet.



---

*Source: 610398-2014-10-27.xml*