Performance Evaluation of Kinetic Codeon Scalar Processors
Takayuki Umeda

Nagoya University, Institute for Space-Earth Environmental ResearchExamples

Eulerian Kinetic (Vlasov) Simulationsfor Space Plasma Studies
Basic equations for collisionless space plasma：• Maxwell equations (for electromagnetic wave propagations)
Computational load
1 ∂E
∇ × B = µ0 J + 2c ∂tless than 0.1%
• Collisionless Boltzmann equation with electromagnetic field(known as Vlasov equation for charged particle motions)f(x, y, z, vx, vy, vz)∂f s∂f s qs∂f s∂t+v⋅∂x+ms(E + v × B ) ⋅∂v=0Operator splitting into three equations [Umeda et al. 2012a]∂f s∂f s+v⋅=0∂t∂x∂f s qs∂f s+E⋅=0∂t ms∂v

6D! ⇒ 5D

∂f s
∂f sqs(v × B ) ⋅=0+
∂t
∂v
ms

（advection in position by v） （advection in velocity by E） （rotation by B）
Non-oscillatory, positive and conservative interpolation is used

System Description[Umeda 2009; Umeda et al. 2012][Umeda et al. 2012b]

405 ~ 4GB
406 ~ 160GB• Hybrid parallelism is adopted toreduce number of processes• Large number of dimensions (up to 6)⇒ Requires huge memory• Length of each loop is short：20-40⇒ number of threads > loop length inmany core environments.
• Multiple loops are thread-parallelized
by loop collapsing of OpenMP[Umeda et al. 2016]

Performance Tuning

1 do n=0,nvz+1
2
do m=0,nvy+1
3
do l=0,nvx+1
4ffi(l)=1.0d0/ff(l,m,n,i,j)
5
end do:30
do l=0,nvx+1
31
gfxy = abs((gfx(l)*ffi(l))* gfy(l))*0.5d032
gfyz = abs((gfy(l)*ffi(l))* gfz(l))*0.5d033
gfzx = abs((gfz(l)*ffi(l))* gfx(l))*0.5d034
gfxyz = abs((gfx(l)*ffi(l))*(gfy(l)*ffi(l))*gfz(l))*inv3
35dfx(l+l0,m,n) = dfx(l+l0,m,n) +(gfx(l)+(gfxyz-gfxy-gfzx)*sx)
36dfx(l+l0,m,n+nv) = dfx(l+l0,m,n+nv)-(gfxyz-gfzx)*sx
37
dfx(l+l0,m+mv,n) = dfx(l+l0,m+mv,n)-(gfxyz-gfxy)*sx
38
dfx(l+l0,m+mv,n+nv) = dfx(l+l0,m+mv,n+nv)+ gfxyz*sx
39
dfy(l,m+m0,n) = dfy(l,m+m0,n) +(gfy(l)+(gfxyz-gfyz-gfxy)*sy)
40
dfy(l+lv,m+m0,n) = dfy(l+lv,m+m0,n)-(gfxyz-gfxy)*sy
41dfy(l,m+m0,n+nv) = dfy(l,m+m0,n+nv)-(gfxyz-gfyz)*sy42
dfy(l+lv,m+m0,n+nv) = dfy(l+lv,m+m0,n+nv)+ gfxyz*sy
43dfz(l,m,n+n0) = dfz(l,m,n+n0) +(gfz(l)+(gfxyz-gfzx-gfyz)*sz)
44dfz(l,m+mv,n+n0) = dfz(l,m+mv,n+n0)-(gfxyz-gfyz)*sz
45dfz(l+lv,m,n+n0) = dfz(l+lv,m,n+n0)-(gfxyz-gfzx)*sz
46dfz(l+lv,m+mv,n+n0) = dfz(l+lv,m+mv,n+n0)+ gfxyz*sz
47
end do48
end do
49 end do

As-is code

There dependences on array access in
the loop iteration at Lines 35—46

30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
5455

• Xeon Phi 7250 – 96GB DDR4Broadwell– Intel Compiler Ver.17.0.1– Option: -ipo -ip -O3 -xMIC-AVX512

KNL• Xeon E5-2697 v4 (dual) – 512GB DDR4
– Intel Compiler Ver.17.0.1– Option: -ipo -ip -O3 -xCORE-AVX2• SPARC64 (Fujitsu FX100) – 32GB HBM– Option: -Kfast,ocl,simd,openmp,parallel• Number of grids: 40*40*40*128*64*2
Fujitsu FX100©Nagoya Univ.
(~28GB)• Elapsed time for 5 time steps is measured

Performance Measurements

3Strong scaling with 𝑁𝑁𝑥𝑥 × 𝑁𝑁𝑦𝑦 × 𝑁𝑁𝑣𝑣𝑥𝑥 × 𝑁𝑁𝑣𝑣𝑦𝑦 × 𝑁𝑁𝑣𝑣𝑣𝑣 = 134 × 70 × 40
Two processes (ion and electron) per compute node

As-is code

Tuned code
Decomposition of most-inner loops ⇒ no loop dependence
Reduction of temporary work arrays

do ll=nvxs-1,nvxe+1
ffi = 1.0d0/ff(ll,mm,nn,ii,jj)gfx(ll) = abs(gfx(ll)*ffi)
gfy(ll) = abs(gfy(ll)*ffi)
gfz(ll) = abs(gfz(ll)*ffi)
end dodo ll=nvxs-1,nvxe+1
dfx(ll+l0,mm,nn) = dfx(ll+l0,mm,nn) + ff(ll,mm,nn,ii,jj)*gfx(ll)*(gfy(ll)*(gfz(ll)*inv3-0.5d0)+(1.0d0-gfz(ll)*0.5d0))*sx
dfy(ll,mm+m0,nn) = dfy(ll,mm+m0,nn) + ff(ll,mm,nn,ii,jj)*gfy(ll)*(gfz(ll)*(gfx(ll)*inv3-0.5d0)+(1.0d0-gfx(ll)*0.5d0))*sy
dfz(ll,mm,nn+n0) = dfz(ll,mm,nn+n0) + ff(ll,mm,nn,ii,jj)*gfz(ll)*(gfx(ll)*(gfy(ll)*inv3-0.5d0)+(1.0d0-gfy(ll)*0.5d0))*sz
end dodo ll=nvxs-1,nvxe+1dfx(ll+l0,mm,nn+nv) = dfx(ll+l0,mm,nn+nv) - ff(ll,mm,nn,ii,jj)*gfx(ll)*gfz(ll)*(gfy(ll)*inv3-0.5d0)*sx
dfy(ll+lv,mm+m0,nn) = dfy(ll+lv,mm+m0,nn) - ff(ll,mm,nn,ii,jj)*gfy(ll)*gfx(ll)*(gfz(ll)*inv3-0.5d0)*sy
dfz(ll,mm+mv,nn+n0) = dfz(ll,mm+mv,nn+n0) - ff(ll,mm,nn,ii,jj)*gfz(ll)*gfy(ll)*(gfx(ll)*inv3-0.5d0)*sz
end dodo ll=nvxs-1,nvxe+1
dfx(ll+l0,mm+mv,nn) = dfx(ll+l0,mm+mv,nn) - ff(ll,mm,nn,ii,jj)*gfx(ll)*gfy(ll)*(gfz(ll)*inv3-0.5d0)*sx
dfy(ll,mm+m0,nn+nv) = dfy(ll,mm+m0,nn+nv) - ff(ll,mm,nn,ii,jj)*gfy(ll)*gfz(ll)*(gfx(ll)*inv3-0.5d0)*sy
dfz(ll+lv,mm,nn+n0) = dfz(ll+lv,mm,nn+n0) - ff(ll,mm,nn,ii,jj)*gfz(ll)*gfx(ll)*(gfy(ll)*inv3-0.5d0)*sz
end dodo ll=nvxs-1,nvxe+1
dfx(ll+l0,mm+mv,nn+nv) = dfx(ll+l0,mm+mv,nn+nv) + ff(ll,mm,nn,ii,jj)*gfx(ll)*gfy(ll)*gfz(ll)*inv3*sx
dfy(ll+lv,mm+m0,nn+nv) = dfy(ll+lv,mm+m0,nn+nv) + ff(ll,mm,nn,ii,jj)*gfx(ll)*gfy(ll)*gfz(ll)*inv3*sydfz(ll+lv,mm+mv,nn+n0) = dfz(ll+lv,mm+mv,nn+n0) + ff(ll,mm,nn,ii,jj)*gfx(ll)*gfy(ll)*gfz(ll)*inv3*sz
end doTuned code

Summary

Scales very well by thread parallelism
 Speedup by performance tuning:–Xeon Phi 7250 (KNL): 2.0 times–Dual Xeon E5-2697 v4 (Broadwell) : 1.5 times–SPARC64 (Fujitsu FX100): 1.15 times
Performance improved by reduction oftemporary work arrays on Broadwell and KNL
 Bottlenecks in scheduling of instructions onSPARX64 XIfx