PIC Simulationin Many-Core Era
Hiroshi Nakashima(ACCMS, Kyoto University)

RI2T ASC Forum © 2014 H. Nakashima

Contents
Supercomputers & Many-Core in Kyoto
 Many-Core (&Exa-Scale) Challenge =Extreme Parallelism (+locality)inter-node: tough intra-node: tougher intra-core: toughest

PIC Simulation on Many-Core
inter-node: OhHelp
 intra-node: particle sorting & cell coloring intra-core: strict particle sorting & ...
Conclusion & Future WorkRI2T ASC Forum © 2014 H. Nakashima

Supercomputers in Kyoto U.
Camphor
XE66300 Abu Dhabi(16 core x 2 socket + 64 GB)
x 940 node
= 30,080 core + 58.75 TB
 300.8 TFlops

Magnolia
XC30

Xeon Haswell(14 core x 2 socket + 64 GB)
x 416 node
= 11,648 core + 26 TB
 428.6 TFlops

Camellia
XC30

Xeon Phi + Xeon((60+10) core + (8+32 GB))x 482 node
= 33,740 core + 18.8 TB
 583.6 TFlops

InfiniBand FDR/QDR

Laurel
GB 8000

Xeon Sandy Bridge
M2090(8 core x 2 socket + 64 GB)
x 601 node (64 w/ GPU)= 9,616 core + 37.56 TB
 242.5 TFlops

Cinnamon
2548X

Xeon Sandy Bridge(8 core x 4 socket + 1.5TB)
x 16 node= 512 core + 24 TB
 10.6 TFlops

InfiniBand FDR

SFA10000

5.0 PB + 3.0 PB
54 GB/sec + 24 GB/sec

Many-Core System in Kyoto U.

RI2T ASC Forum © 2014 H. Nakashima

Cray XC30 Camellia
GDDR58GBBlade

opt:5x3x12.5Gbps BP:1x15x14Gbps

chassisDDR3-1600
32GB

cabinet#0cabinet#1

copper:5x3x14Gbps Xeon E5-2670v2

Xeon Phi 5120D
1.053GHz x 16 x 60
= 1011GFlops (60+10)C x 482 = 33740C
 583.6TFlops
 18.8TB

 4.6TB/s

2.5GHz x 8 x 10
= 200GFlops

cabinet#2copper:48x5.25GB/s
opt:8x4.68GB/s

RI2T ASC Forum © 2014 H. Nakashima

Many-Core & Exa-Scale Challenge

Even with 0.5PF+ Camelliahuge parallelism of482node x 120thread x 16DPFlop = 925,440DPFlop

A 30PF/1000PF system could have
10,000 x 100-200
 100,000 x 150-300
x 32x 64= 32-64M= 1000-2000M

Which is toughest ininter-node / intra-node / intra-core ?
We have to play with MPI for inter-node.
 OpenMP provides us of easy parallelization.
 We can rely on compilers for vectorization.
Really?

Many-Core & Exa-Scale Challenge

RI2T ASC Forum © 2014 H. Nakashima

Inter-Node is Toughest?
It's tough but not toughest because;We've written (scalable) MPI codes.
 We've learned scalability issues.


multi-dimensional/paradigmatic decomposition
 minimizing frequency/amount of collectives
 other bottlenecks e.g., initialization, I/O, etc.




At least we have forerunners on (e.g.) K.

c.f. Latency hiding by comp/commoverlapping is much tougher, but;Your code cannot be slower than 2x without it(or cannot be faster than 2x with it).
 Unless we really need latency-hiding algorithms(not implementations).
Many-Core & Exa-Scale Challenge

RI2T ASC Forum © 2014 H. Nakashima

Intra-Node is Easier?
It was easier than inter-node in goodold 101-scale multithreading.
Simple 1D decomposition with"#pragma omp (parallel) for" was OK.
1
 10 -order overhead was acceptable.


Manycore x HT changed the story.
2403 subproblem is too small/large in terms of
performance/memory-size.
 MPI-like coding is not only cumbersome but alsoeasy to bring performance bugs;


because we can easily do a stupid thing which we
couldn't do with MPI.

Best decomposition for HT may be different frombest for ordinary multithreading.

Many-Core & Exa-Scale Challenge

RI2T ASC Forum © 2014 H. Nakashima

Intra-Core is Toughest?


Definitely YES because;
SIMD ≠ vector due to lack of list-vectors, large
memory bandwidth, etc.
 We have to rely on compilers.


OpenMP made us free from auto-parallel compilers.

Therefore we'll have to do ...
Rewrite innermost kernels applying loop-split,loop-interchange, loop-xxx by hand taking care ofaccess locality.
 Sometimes redesign data structure to conformwith 1D-innermost parallelism.
(e.g. array-of-structure to structure-of-array)

RI2T ASC Forum © 2014 H. Nakashima

PIC Simulation on Many-Core

What's that?
Simulate motion of a huge
number of (e.g., 1011) ion/
electron particles in alarge scale (e.g., (103)3)
discretized electro-magnetic field.
 particle push: A particle is accelerated by Lorentz
force calculated from discretized E&B around thecell (voxel) where the particle resides, current scatter: while its motion contributes to thecurrent around the cell, field solve: which then affects E&B by Maxwell'slow.
RI2T ASC Forum © 2014 H. Nakashima

PIC Simulation
More formally ...

~ ≡ [0, N ]3 , ∆x ≡ X /( N + 1)
Ω ≡ [0, X ]3 , Ω~ | x / ∆x − x~∈Ω
~ < 1}
δ x ≡ {x∞∀( x i , vi , qi ) ∈ Π : vi += f ( E (δ x i ), B(δ x i ), x i , vi , qi ); x i += vi ;~ : J(x~∈Ω~ ) = ∑ g( x − v , v , q ) − ∑ g( x ,− v , q );∀xxi − vi ∈δ−1 ~xiiiixi ∈δ−1 ~xiii~ : E( x~∈Ω~ ) += (∇ × B( x~) / µ − J ( x~ )) / ε ; B( x~ ) −= ∇ × E ( x~ );∀x~Ω, ΩE (δ xi ), B(δ xi )J (δ ( xi − vi ))J (δ xi )RI2T ASC Forum © 2014 H. Nakashima

PIC Simulation
More procedurally ...
for(t){
for(k)for(j)for(i) J[k][j][i]=0;
for(p){
v[p]=lorentz(...); x[p]+=v[p];scatter(x[p]-v[p],v[p], q[p],J);
scatter(x[p],-v[p],-q[p],J);}
filed_solve(E,B,J);}
void scatter(x,v,q,J) { (i,j,k)=x;J[k ][j ][i ]+=...; J[k ][j ][i+1]+=...;...
J[k+1][j+1][i ]+=...; J[k+1][j+1][i+1]+=...;}

PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Inter-Node: OhHelp

What's that?

03132333av.# of particles

02 03 12 21 22 23 3200 3132103301
11 20 21 01 31 02
1130
132200102030

Dynamic load balancer forforparticle simulations tocope with extreme
imbalance of particle density.
 A process is responsible for its own primarysubdomain and a secondary subdomain of other
having too many particles.
x2865
16000
 Good scalability up to
uniform
12000
x2720
4096 processes (andcongestedmore) even for the case 8000
in which particles
4000
congest in a small
0region.
33 00 32 01 30 10 13 03 23 20 31 02 11 21 12 22

av.# of particles

33 00 32 01 30 10 13 03 23 20 31 02 11 21 12 220

1000

2000

#processes

3000

4000

PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Node: Was Easy


OpenMP code for 101-scale threadingfor(t){#pragma omp parallelme=omp_thread_num();for(k)for(j)for(i) JJ[me][k][j][i]=0;#pragma omp for
for(p) {v[p]=lorentz(...); x[p]+=v[p];scatter(x[p]-v[p],v[p], q[p],JJ[me]);scatter(x[p],-v[p],-q[p],JJ[me]);}#pragma omp forfor(k)for(m)for(j)for(i)J[k][j][i]+=JJ[m][k][j][i];...

PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Node: Is Tough with KNC


Temporary array JJ is too large.
|JJ|=N3 x 3 x #threads3
 |(x,v)|=N x 6 x average-particle-density
 #threads=240, average-particle-density=100|JJ|=N3 x 720 > |(x,v)|=N3 x 600
 Making memory space >2x is prohibitive for KNC'ssmall memory.
Scanning large JJ twice in a time step;
Makes required bandwidth >3x.
 Limits performance due to memory bandwidthbottleneck even with GDDR5.
PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Node: Particle SortingParticle sorting according to position;
Improves access locality of E, B & J.
 Enables us to hold all particles in a set of cells.


Sorting can be non-strict.
May ignore some of coordinates (e.g., sort only byz-coordinate).
 May assume coordinate values in a range equal(e.g., regarding all particles in a cell-block equallypositioned).
 May be done occasionally (e.g., every t time steps).


Sorting can be O(N) because we know
the range of values.

PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Node: Cell-Block ColoringCell-block coloring with 2k colors for
k-D decomposition for multithreading;Makes current-scatter in blocks of same color
conflict-free, if block-size≥3 and we hold allparticles in each block.

Sorting can be non-strict, e.g.;
every t steps if block-size≥3+t.
 if "omp atomic" is fast for non-conflict write.
PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Node: KNC PerformancePreliminary experiment setting
963 x 100 particles/cell, uniformly distributed
 2D-decomposition, 4 blocks/threadPerformance
Looks good.
 But not actually.
HT=180200

x57.3

106 particles/sec

106 particles/sec



HT=260

HT=3
HT=440200150100
500020#cores4060
HT
GFlops

KNC
21011

Haswell
2515

Ivy Bridge
2200

PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Core: Why So Slow?
Innermost loops is too complicated.
Read/write indices of E/B/J in an iteration looks
random for compiler.
 Tough to SIMDize 8x3x2 reads of E&B.
 Impossible to SIMDize 8x3x2 reads/writes of J.
E[][][]/B[][][]
for(p in a block)48

48v[p]=lorentz()
v[p]=lorentz()
x[p]+=v[p]x[p]+=v[p]scatter(x[p]-v[p],…) scatter(x[p]-v[p],…)scatter(x[p],…)scatter(x[p],…)48

48

J[][][]

48
v[p]=lorentz()x[p]+=v[p]scatter(x[p]-v[p],…)scatter(x[p],…)48

PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Core: Make Loop SIMDizable

Sort particles cell-wise & strictly.
Indices of E, B & J (before push) are common for
all particles in a cell.
 Therefore elements accessed in the loop for a cell
can be placed in a set of 72(!!) scalar variables.
Split particle loop into 3 loops;
To update v[p] referring to scalarized E&B.
 To scatter currents into scalarized J and thenupdate x[p], for reducing the number of commonscalars.
 To scatter currents after sorting, so that J's
indices are common for particles in their newresidence cell and thus J is scalarized.
PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Core: Make Loop SIMDizablefor(each block & cell)
E[][][]/B[][][]
48

scalarized E/Bfor(p in a cell)v[p]=lorentz()
for(p in a cell)
scatter(x[p],…)x[p]+=v[p]v[p]=lorentz()

v[p]=lorentz()

scatter(x[p],…)x[p]+=v[p]scatter(x[p],…)x[p]+=v[p]24scalarized J

J[][][]

sort()
for(each block & cell)for(p in a cell)
scatter(x[p],…)scatter(x[p],…)scatter(x[p],…)24

scalarized J

J[][][]

PIC Simulation on Many-Core

RI2T ASC Forum © 2014 H. Nakashima

Intra-Core: KNC Performance
Large improvement & good speedup.
 And actually good now.
600
HT=1
HT=2500

500

HT=3400

106 particles/sec

106 particles/sec

600
x6.9

x51.5

HT=4old300
200

400
300
200

100

100

00020

#cores4060

HT
GFlops

KNC
31011Haswell2515

Ivy Bridge
2200

RI2T ASC Forum © 2014 H. Nakashima

Conclusion & Future Work

What we've learnt.
Simple porting or simple large-scale multithreading may be far from sufficient to have good
performance on KNC (as anticipated).
 SIMD-aware kernel recoding/redesign may benecessary in KNC playground (as anticipated).
What we'll have to do.
Design efficient, locality- & SIMD-aware, sorting.
 Make two-phase current scatter more localityconforming (e.g., by tiling).
 More generally, improve our DSL framework oflocal-view coding scheme to generate the tuned
code (semi-)automatically.
