Convergence property of Krylov subspace methods forasymmetric linear system on fusion plasma simulation code
GT5DSusumu Yamada1,3, Yasuhiro Idomura1, Toshiyuki Imamura2,3, Masahiko Machida1,312

Japan Atomic energy Agency

The University of Electro-Communications
3CREST (JST)

1. Introduction
The gyrokinetic toroidal five-dimensional Eulerian code
GT5D [1] is one of the first-principle simulation codes for thefusion plasma turbulence. The code has been parallelized andported on some massively parallel platforms, and then, it has
been confirmed that the code achieves good parallel efficiency
on them. In GT5D code, the second-order additive
semi-implicit Runge-Kutta (ASIRK) method [2] is adopted forthe time-evolution simulation. Therefore, we need to solve thelinear equations, whose coefficient matrices are asymmetric. In
the original code, the generalized conjugate residual (GCR)
method without preconditioning has been so far utilized. Thus,
it is expected that the convergence property can be improved bysubstituting a more suitable preconditioned solver. In thisresearch, we examine some preconditioned Krylov subspacemethods and measure their convergence properties. Moreover,
we evaluate the parallel performance of the solvers.
2.

Convergence Property of Krylov Subspace Methods

2.1 Preconditioners
Here, we discuss the effect of preconditioning for the linear
equations. At first, we examine the convergence property of theincomplete LU (ILU) decomposition. The ILU decomposition
can improve the convergence property for various linear
equations. However, when we employ the ILU decomposition,
the convergence property in GT5D code becomes worse. Thereason is that the ILU decomposition does not approximate thecoefficient matrix well, since the diagonal elements of thecoefficient matrix are small as compared with the off-diagonalelements.
Next, we propose the preconditioning based on the SOR
iteration method. The preconditioning is executed according tothe following formulaxin 1  xin  i 1bi   mij x jmii j 1n 1N  mij x nj  ,j iwhere mij is the ( i, j)-element of preconditioner M, which is anapproximate matrix of the coefficient matrix. Table 1 shows theeffect of the preconditioning based on the SOR iteration. The
problem size is 60×60×32×28 in xyzv-space. Here, since thecoefficient matrix of the linear equation is given by 17 stencil
grids in 4-dimensional space, the preconditioner can becomposed of the elements in each direction. Therefore, we
compare the convergence property for the some preconditioners
covering all or partial dimensions. The result shows that thepreconditioner covering all directions is the best and the effectof the v-direction does not almost exist.
Table 1. Relationship between the dimensions covered by thepreconditioner M and the number of iterations forpreconditioned GCR method.
Preconditioner
ω
Number of(covering dimensions)
iterations
None－
466
x, y-direction－
Not converge
x, z –direction
0.047
289
y, z –direction
0.043
282
z, v –direction
0.049
347
x –direction－
Not converge
z –direction
0.049
354
v –direction
1.8
466
x, y, z, v –direction
0.054
181
(all dimensions)
2.2 Convergence Property of Preconditioned Krylov
Subspace Methods
Here, we examine the effects of the preconditioners covering
all or some dimendions for GCR method, Bi-CGstab method,
and GMRES method.
GCR method
We show the convergence property in Figure 1. The result
demonstrates that the number of iterations decreases from 466
to 181 by utilizing the preconditioning.
Bi-CGstab method
Figure 2 shows the convergence property of the Bi-CGstabmethod. The convergence of the Bi-CGstab method without
preconditioning is very slow (the method can converge by
about 5000 iterations). The preconditioning yields a large
convergence accelaration.
GMRES method
Figure 3 shows the convergence property of the Bi-CGstabmethod. Here, the restart parameter is 4. The number of
iterations decreases to about 1/4 by utilizing thepreconditioning.
These results demonstrate that the preconditioning based on
the SOR iteration method is suitable for this equation.

Number of iterations
0200

400

600

800

1000

1.00E+00
1.00E-01

GCR

1.00E-02
1.00E-03

PGCR

1.00E-04

Residual

1.00E-05
1.00E-06
1.00E-07
1.00E-08

3.

Numerical Experiments
We examine the parallel performance of the preconditioned
solvers on 256 cores of Fujitsu BX900 in Japan Atomic Energy
Agency. The problem size is 160×160×32×128 in xyzv-space.
Since table 1 shows that the effect of the preconditioning forthe ν-direction does not almost exist, we employ thepreconditioners covering the all directions and x,y,z directions.
We show the number of the iterations and the elapsed times inTable 2. The result indicates that the effects of bothpreconditioners are almost the same. Since calculation cost of
the x,y,z-direction preconditioner is less than that of
all-direction one, the x,y,z-direction preconditioner can be
performed faster.

1.00E-09
1.00E-10
1.00E-11
1.00E-12
1.00E-13
1.00E-14
1.00E-15
1.00E-16

Figure 1. Convergence property of GCR method.

Table 2. Number of iterations and Elapsed time on 256 cores ofBX900. P and Pxyz means all-direction preconditioner andx,y,z-direction one, respectively. Since Bi-CGstab method
without preconditioning does not converge within 1000
iterations, its result is omitted.
methodω

Number ofElapsed timeiterations

(sec)

641
260
256
144
141
210
60
60

110.7195107.2379
104.7527
110.8538
97.9747
208.4488
130.5454
118.5158

Number of iterations
0200

400

600

800

1.00E+01

1000

GCRPGCR
PxyzGCR
PBi-CGSTAB
PxyzBi-CGSTABGMRESPGMRES
PxyzGMRES

1.00E+00
1.00E-01
1.00E-02

Residual

1.00E-03
1.00E-04
1.00E-05
1.00E-06
1.00E-07

Bi-CGSTAB

1.00E-08
1.00E-09

PBi-CGSTAB

1.00E-10

1.00E-11－
0.037
0.040
0.038
0.038
－0.0470.047

1.00E-12
1.00E-13

4.

1.00E-14
1.00E-15

Figure 2. Convergence property of Bi-CGstab method. The
convergence of the method without preconditioning is very
slow. On the other hand, the preconditioned method can solve
the equation very fast.
Number of iterations
050

100

1502001.00E+00
1.00E-01
1.00E-02

1.00E-03

Residual

1.00E-04
1.00E-05

GMRES

PGMRES

1.00E-06
1.00E-07

1.00E-08
1.00E-09
1.00E-10
1.00E-11
1.00E-12
1.00E-13
1.00E-14
1.00E-15

Figure 3. Convergence property of GMRES method.
The restart parameter is 4. The number of iterations decreases
to about 1/4 by utilizing the preconditioning.

Conclusion
We examined the convergence property of some Krylov
subspace methods for the linear equation on GT5D code. Theresult showed that the Bi-CGstab method with thepreconditioning based on the SOR iteration is the fastest solver.
Acknowledgement
This work is supported by the MEXT grant for HPCI
Strategic Program Field No.4: Next-Generation Industrial
Innovations. This research was partially supported by CREST,
JST and a Grant-in-Aid for Scientific Research from JSPS
(Grant No. 23500056).
Reference[1] Y. Idomura, M. Ida, T. Kano, N. Aiba, and S. Tokuda.
Conservative global gyrokinetic toroidal full-five-dimensionalvlasov simulation. Computer Physics Communications, Vol.
179, pp. 391-403, 2008.
[2]
Xiaolin Zhong. Additive semi-implicit runge-kutta
methods for computing high-speed nonequilibrium reactive
flows. Journal of Computational Physics, Vol. 128, pp. 19-31,
2007.