Fourth Workshop on Accelerator Programming Using Directives (WACCPD),
Nov. 13, 2017

Implicit Low-Order Unstructured Finite-Element
Multiple Simulation Enhanced by Dense Computation
using OpenACC
Takuma Yamaguchi, Kohei Fujita, Tsuyoshi Ichimura,
Muneo Hori, Lalith Maddegedara, Kengo Nakajima1Introduction• Contribution of HPC to earthquake mitigation highly anticipated from society• We are developing comprehensive earthquake simulation that simulate all phases ofearthquake disaster by use of full K computer system• Simulate all phases of earthquake by speeding up core solver• SC14 Gordon Bell Prize Finalist, SC15 Gordon Bell Prize Finalist & SC16 Best Poster & SC17Best Poster Finalist• Ported this solver to GPU environment using OpenACC in WACCPD 2016 (Best Paper)• Today’s topic is enhancement of this GPU solver, and report performance on Pascal andVolta GPUs

Earthquake disaster process

K computer: 8 core CPU x 82944 node systemwith peak performance of 10.6 PFLOPS

2Comprehensive earthquake simulationc) Resident evacuationa) Earthquake wave propagation

0 km
-7 km

Two million agents evacuating to nearest safe siteIkebukuro

Earthquake

Post earthquake

UenoShinjuku
Tokyo station

ShinbashiShibuyab) City response simulation

Large finite-element simulation enabled
by developed solver

3

City simulation

Visualized by CYBERNET SYSTEMS CO., LTD4

Target problem: Earth’s crust deformation problem• Compute elastic response to given fault slip
• Many case analysis required for inverse analyses and Monte
Carlo simulations• Compute using finite-element method: solve large
matrix equation many times0• Involves many random data access & communication• Difficulty of problem• Attaining load balance & peak-performance & convergency ofiterative solver & short time-to-solution at same time
• Smart use of compute precision space, constraints in solver
search space according to physical solution space required

Ku = fEarth’s crust deformation
problem

Outer force vector

Unknown vector with up to 1 trillion degrees of freedom

Sparse, symmetric positive definite matrix

5Designing scalable & fast finite-element
solver• Design algorithm that can obtain equal granularity at O(million) cores• Matrix-free matrix-vector product (Element-by-Element method) is promising:Good load balance when elements per core is equal• Also high-peak performance as it is on-cache computation• Combine Element-by-Element method with multi-grid, mixed
precision arithmetic, and adaptive conjugate gradient method• Scalability & peak-performance good (core computation kernels are Elementby-Element), convergency good, time-to-solution good
Element-by-Element method+=

Element #0+=Element #1

f = Σe Pe Ke PeT u[Ke is generated on-the-fly]

Keu

…

…

f

Element #N-1

6

Solver algorithm
Solving preconditioning matrix (single precision)
Inner loop level 2
CG loop

Solve system roughly using CG solver

Second ordered
tetrahedron

Computations ofouter loop

Use x1←P12 x2 as
initial solution

Solve system roughly using CG solver

Outer loop
Inner loop level 0

Use x0←P01 x1 as
initial solution

Linear tetrahedron

Geometriccoarsening

Equation to be solved(double precision)

Inner loop level 1

Algebraiccoarsening

Solvingpreconditioning
matrixCoarsened equation(P12TA1P12)x2=P12Tb1

Solve system roughly using CG solver

K. Fujita, T. Ichimura, K. Koyama, H. Inoue,
M. Hori, L. Maddegedara, Proceedings of
the Platform for Advanced Scientific
Computing Conference (PASC), June 2017

Use for preconditioner of outer loop

Second orderedtetrahedron
7Performance on K computer• Developed solver significantly faster than

• PCGE (standard CG solver algorithm; preconditioning with 3x3 block diagonal matrix)• SC14 Gordon Bell Prize finalist solver (base solver for WACCPD 2016 GPU solver)• Use this as a base for GPU solver
1800

9,216

18,43236,86473,728 147,456 294,912# of CPU cores

1400

200

146.5 s (20.6%)
319.7 s (20.4%)
1629 s (17.6%)400148.7 s (20.1%)
308.5 s (20.5%)1573 s (17.9%)600

151.3 s (20.1%)
311.8 s (21.3%)
1680 s (18.3%)

800

144.5 s (20.9%)
292.7 s (21.4%)
1572 s (18.1%)

1000

140.0 s (21.2%)
277.5 s (21.6%)
1546 s (18.1%)

1200
138.9 s (21.3%)
272.0 s (21.9%)
1375 s (18.5%)

elapsed time of solver [s]1600

Floating-point
arithmetic
efficiency to peakElapsed time

• 94.8% scalability from9216 to 294912 cores• 4 times peak performanceof HPCG benchmark(HPCG on K computer:
5.3% in double precision)

08

Introduction of GPU computations• Further speedup of the simulation by introducing GPUs
• Good load balance, Reduced computation cost & data transfer size isalso beneficial for GPUs
• High performance can be obtained using OpenACC
with low development cost

• GPU architecture is different from CPU architecture• Latency bound especially when we conduct random memory access• Relatively smaller cache size

• Finite-element applications tend to be memory bandwidth bound➔Simple porting of the CPU code is not sufficient9Key kernel: Element-by-Element kernel• Most costly kernel; involves data recurrence• Algorithm for avoiding data recurrence on CPUs
• Use temporary buffers per core & per SIMD lane
• Suitable for small core counts with large cache capacity• Algorithm for avoiding data recurrence on GPUs
• Last year, we developed algorithm using atomics and achieved high performance• However, random access becomes bottleneck…

Element-by-Element method+=

Element #0+=Element #1

Data recurrence(add into same node)

Keu…

…

f

Element #N-1

10

Performance in simple porting
DOF: 125,177,217, # of elements: 30,720,000

Outer

Inner level 0

Inner level 1

Inner level 2

Computational Environment
K computer

Reedbush-H# of nodes

2010CPU/node

1x
SPARC64
VIIIfx2 x Intel
Xeon E52695 v4--

2 x NVIDIA
P100

GPU/node

Hardware128
peak FLOPSGFLOPS/process

5.30
TFLOPS

Memorybandwidth/process

732 GB/s

64 GB/s

39.00

K computer
DP_FLOPS efficiency (%)
MEM. Efficiency (%)

Reedbush-H

24.91
22.81

21.23
16.78

7.75
010

9.54
43.81

21.47
23.18

1/5.0
2030

40

Elapsed time of solver [s]• Simple porting achieved 5.0 times speedup• However, there is some room for improvement• Memory bandwidth is 11 times larger
11Strategy for Introduction of OpenACC• To attain optimal performance, algorithm/implementation
suitable for GPUs should differ from that for CPUs
Thereby, we
1. Modify the solver algorithm to suit the GPU architectures
2. Port the solver to GPUs using OpenACC

12

Modification of Algorithm for GPUs
• Reduce random memory accesses• Target applications (Inverse analyses, Monte Carlo method etc.)
solve many systems of equations• Same stiffness matrix• Different right-hand side input vectors

Element-by-Element kernel
Element #0• Multiple equations at the same time
K u1 , u2 , u3 , … , u16 𝑇 = f1 , f2 , f3 , … , f16 𝑇
Instead of Ku1 = f1 , Ku2 = f2 , Ku3 = f3 , …

k…
f1 f2 f3 f4…

e

Element #1…

ke…

u1 u 2 u3 u413

Introduction of OpenACC – 1/3
Control of data transfer
CPU computation

Read data
Input vector f!$acc data copy(u, …) copyin(f, …) GPU computation
Solve u=K−1 fComputation (i.e, EBE)!$acc update host(err)Check convergence!$acc end data

Output vector u

Data transfer between CPU and GPUs
is minimized in the solver• Only in convergence check part
• GPU Direct is used for
MPI communication
14

Introduction of OpenACC – 2/3
Insertion of some directives forparallel computation!$acc parallel loop collapse(2) present(…
do i_ele = 1, n_element
do i_vec = 1, 16
cny1 = connect(1, i_ele)…
cny10 = connect(10, i_ele)
u0101 = u(i_vec, 1, cny1)…
u1003 = u(i_vec, 3, cny10)
Ku01 = …

…

Example for Element-by-Element
multiplication• Assign 16 threads for one
element• Introduce atomic functions toavoid data race

Ku30 = …
!$acc atomicr(i_vec, 1, cny1) = r(i_vec, 1, cny1) + Ku01

…

12
3
4
5678910
1112
13
14
15
161718
1920
21
2223!$acc atomicr(i_vec, 3, cny10) = r(i_vec, 3, cny10) + Ku30
enddo
enddo
15!$acc end parallel

Introduction of OpenACC – 3/3
Minor tuning for OpenACC parameters• The allocation of gang, worker and vector• The length of vector
Optimize fine-grain control of parallelism(Not large effect on performance)

16Performance of the proposed solver
DOF: 125,177,217, # of elements: 30,720,000

Outer

Inner level 0

Inner level 1

Inner level 2

Computational Environment
K computer

Reedbush-H# of nodes

2010

CPU/node

1x
SPARC64
VIIIfx2 x Intel
Xeon E52695 v4--

2 x NVIDIA
P100

GPU/nodeHardware
128
peak FLOPSGFLOPS/process

5.30
TFLOPS

Memorybandwidth/process

732 GB/s

64 GB/s

K computer
1 vector

39.00
DP_FLOPS efficiency (%)
MEM. Efficiency (%)

Reedbush-H
1 vector(24.9122.81

21.23
16.78

7.75

9.54
43.81

21.47
23.18

1/5.0

K computer
computerK
16 vectors
vectors1626.43)

1/9.6

Reedbush-H
16 vectors

2.75

1/14.20102030

Elapsed time of solver, per vector [s]1740The speedup of each kernel
Elapsed time per vector(s)Kernel1 vector

16 vectors

Speedup

1st order (FP32)

0.948

0.584

1.62

2nd order (FP32)

0.687

0.401

1.71

2nd order (FP64)

0.044

0.025

1.78

SpMV

1.465

0.091

16.10

Dot product

0.213

0.522

0.41

Total

7.75

2.75

2.82

Element-byElementcomputation• Reduction in random memory access in EBE kernels• Total computation time for SpMV is constant• Bound by reading global matrix for memory• Dot product is not efficiently computed• OpenACC cannot use arrays for reduction option• Using scalars (tmp1,tmp2,…,tmp16) causes stride memory access18

Weak Scaling
Reedbush-H: P100 GPU x 240
ModelDOF# of elements

PCGE
iterations# ofGPUs

Outer

Inner level 0

Inner level 1

Inner level 2

12 + 311 + 1300 + 2161

No. 1

125,177,217

30,720,000

4,928

20 No.1

1.57

16.76

21.2147.47

7.93

12 + 311 + 1346 + 2477

No. 2249,640,977

61,440,0004,943

40 No.2

1.60

16.95

22.05

49.63

9.03

11 + 280 + 1312 + 1924

No. 3496,736,817

122,880,0004,901

80 No.3

1.50

15.54

21.61

6.99

45.64

11 + 281 + 1430 + 2367

No. 4992,038,737

245,760,0004,905

160 No.4

1.51

15.83

23.68

8.97

49.99

23.27

8.95

49.89

11 + 280 + 1341 + 2201

No. 5 1,484,953,857

368,640,0004,877

240 No.5

1.540

16.13

10203040

Elapsed time of solver [s]5019Performance in using V100 GPUs
DOF: 38,617,017, # of elements: 9,440,240
Outer

Computational Environment
Reedbush-H

DGX-1 (P100/V100)# of nodes

41CPU/node

2 x Intel Xeon
E5-2695 v42x
Intel Xeon E5-2698 v4

2 x NVIDIA
P100

8 x NVIDIA 8 x NVIDIA
P100
V100

GPU/node

Hardwarepeak
FLOPS/ process

5.30
TFLOPS

Memorybandwidth/process

732 GB/s

5.30
TFLOPS

7.5
TFLOPS

Inner level 0

Inner level 1

Inner level 2

Reedbush-H(P100 GPU)

19.21

DGX-1(P100 GPU)

17.28

DGX-1(V100 GPU)

6.86
05101520Elapsed time of solver [s]732 GB/s

900 GB/s

Higher performance than expected fromhardware capability• improved L1/L2 caches
20Application Example
Estimation of coseismic slip distribution in 2011 Tohoku Earthquake
DOF: 409,649,580# of input vectors: 368 = 23 sets × 16 vectors
Computation Environment: 64 x P100 GPUs (32 nodes of Reedbush-H)
Computation time: 828s for 23 sets of analyses
Target Domain(29 times better in performance than previous studies)(i) Whole FE modelx(ii) Close-up viewy1192 km

792 km
400 kmz(iii) Material properties
Crust layer
Upper-mantle layer
Philippine Plate
Pacific Plate

Vp (m/s) Vs (m/s) ρ (kg/m3)
5664
3300
2670
8270
4535
3320
6686
3818
2600
6686
3818
2600

21

Application Example
Estimation of coseismic slip distribution in 2011 Tohoku Earthquake
DOF: 409,649,580# of input vectors: 368 = 23 sets × 16 vectors
Computation Environment: 64 x P100 GPUs (32 nodes of Reedbush-H)
Computation time: 828s for 23 sets of analyses
Estimated Distribution(29 times better in performance than previous studies)(i) Whole FE modelx(ii) Close-up viewy1192 km

792 km
400 kmz(iii) Material properties
Crust layer
Upper-mantle layer
Philippine Plate
Pacific Plate

Vp (m/s) Vs (m/s) ρ (kg/m3)
5664
3300
2670
8270
4535
3320
6686
3818
2600
6686
3818
2600

22

Conclusion• Accelerate the unstructured low-order finite element solvers byOpenACC• Design the solver appropriate for GPU computations• Port the codes to GPUs

• Obtain high performance with low development costs• 14.2 times speedup on P100 GPUs from the original solver
on CPU-based K computer• Computation in low power consumption• Improvement in reliability of earthquake simulation• Many-case simulation within short time23