Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2020

jh200051-NAHI

Scalable Multigrid Poisson solver for AMR-based CFD applications
in Nuclear Engineering
Naoyuki Onodera（Japan Atomic Energy Agency）
We develop a multigrid preconditioned conjugate gradient (MG-CG) solver for thepressure Poisson equation in a two-phase flow CFD code JUPITER. The MG-CG
solver is redesigned to realize efficient CFD simulations including complex
boundaries and objects using an adaptive mesh refinement (AMR) method based ona block-structured Cartesian grid system. Here, we propose a new MG preconditioner
with the cache reuse successive over relaxation (CR-SOR) smoother, which has higharithmetic intensity using the shared memory and enables continuous memory accesson the block-structured grid. We measured the kernel performance of the MG-CG
solver on GPU and CPU on the TSUBAME supercomputer, and A64FX on the FLOW
supercomputer. The performances on the GPU and CPU are reasonable, but there are
significant performance degradation on A64FX. The numerical experiments are
conducted for two-phase flows in a fuel bundle of a nuclear reactor. Thanks to theblock-structured data format, grids inside fuel pins are removed without performancedegradation, and the total number of grids is reduced to 2.26×109, which is about 70%of the original uniform Cartesian grid system. The MG-CG solver with the CR-SOR
smoother reduces the number of iterations to less than 9% of the original
preconditioned CG method, leading to 5.9-times speedup on the TSUBAME
supercomputer. In the strong scaling test, the MG-CG solver with the CR-SOR
smoother is accelerated by 2.1 times between 64 and 256 GPUs.

1.

Basic Information

・ Yuichi Asahi (JAEA, Porting codes on CPU/GPU
platforms)(1) Collaborating JHPCN Centers
・ The University of Tokyo・ Takashi Shimokawabe (Tokyo, Advice for largescale computations on Oakbridge-CX)・ Tokyo Institute of Technology

・ Satoshi Oshima (Nagoya, Advice for large-scale

・ Nagoya University

computations on Flow)(2) Research Areasn

Very large-scale numerical computation

2.

(3) Roles of Project Members
・ Naoyuki

Onodera(JAEA,(1) Research purpose

Representative,

This project aims at developing exascale computing

Development of a Poisson solver)
・ Audit

Edouard(CEA

France,technologies for multi-scale CFD simulations on

Deputy-exascale platforms. For this purpose, we promote a

Representative, Advice and support)
・ Takayuki

Aoki(TokyoTech,

Purpose and Significance of Researchcollaborative research with respect to the following

Deputy-three main subjects, a) High performance sparse matrix

Representative, Advice and support)

solvers on accelerated computing platforms, b)

・ Yasuhiro Idomura (JAEA, Advice and support for

Adaptive Mesh Refinement (AMR) methods for multi-

large-scale computation on TSUBAME)

scale turbulence problems with complicated boundaries,

・ Yuta Hasegawa (JAEA, Development of a Poissonc) Performance portability of CFD codes based on

solver)advanced frameworks and programming models.

・ Susumu Yamashita (JAEA, Development of a

The target CFD applications in this project are thenuclear engineering code)1

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2020

3D multi-phase multi-component thermal hydraulic
CFD code JUPITER [Yamashita et al., Nuclear
Engineering and Design 2017] and the 3D compressiblehydrodynamic astrophysics CFD code ARK [Padioleau
et al., Astrophysical Journal 2019]. JAEA has promoted
the development of high-performance pressure Poissonsolvers based on various Communication-Avoiding
Fig.1: Strong scaling tests up to 7,680 GPUs on theSummit. A large JUPITER matrix with N = 1,280 × 1,280
× 4,608 is computed using the Preconditioned CG method
and the Preconditioned Chebyshev Basis CommunicationAvoiding CG method.

(CA) algorithms on the Oakforest-PACS [Idomura,
ScalA18@SC18]

and

on

the

Summit[Ali,

ScalA19@SC19] (Fig.1), and their strong scaling wasdemonstrated up to ~8,000 KNLs/V100s. In the
FY2019 project, the development of a new AMR
framework for JUPITER was conducted. On the other

JUPITER. If such an efficient implementation works forhand, CEA has developed ARK based on an AMR-

AMR-CFD simulations, it has a great impact on various

based compressible fluid model, and its performanceengineering and science fields.

portability between KNL and V100 was demonstratedb) AMR methods for multi-scale CFD problems

using a programing model with MPI+Kokkos. ARK

AMR methods are promising solutions for multi-scalerequires high performance matrix solvers for gravity

turbulence problems. Their parallel implementations

Poisson solvers and implicit matrix solvers.

have been designed to keep the balance of computing

In this project, we accelerate the development ofloads, e.g., via space filling curves. However, oncommon exascale computing technologies for AMR-

accelerated computing platforms, the communication

CFD simulations by sharing the knowledges on highcost often becomes comparable or exceeds theperformance matrix solvers, performance portability,

computing cost. In such a situation, one needs to design

and AMR methods. The developed technologies willa new parallel implementation by taking account of bothcontribute to various science and engineering CFD

computation and communication costs. In this work, weapplications.

develop a new performance model of AMR-CFD(2) Significance of research

simulation, which is of critical importance for variousa) High performance sparse matrix solvers on

exascale AMR applications.

accelerated computing platformsc) Performance portability on multi-platforms

Sparse matrix solvers for the pressure Poisson equation,

Performance portability between CPU and GPU has

the gravity Poisson equation, and implicit time

been a long-standing issue. In addition, several new

integration are critical bottlenecks for high performance

CPU and GPU architectures are emerging towards

CFD simulations. In this project, we extend our CA

exascale supercomputers. This makes the performancealgorithms such as CA multi-grid (CA-MG) solvers andportability issue more serious. In this work, first, we

CA-Krylov solvers to block (or patch) structured AMR

port miniapps from JUPITER and ARM via native

grid data, and apply them to various nuclear engineering

programing models such as CUDA and OpenMP andproblems and astrophysics problems in JUPITER and

via more portable programing models such as Kokkos,

ARK, respectively. By analyzing required precision in

and investigate performance bottlenecks in the lattereach problem, we explore optimum mixed precision

approach. From this performance study, we clarify theapproaches, which were already implemented inmost efficient porting strategy and apply it to JUPITER
2

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2020

and ARK.

accurate centered finite difference scheme with sevenstencils. Since the matrix is symmetric and diagonal-

3.

Significance as JHPCN Joint Research Project

dominant, the Poisson solver is normally computed

The goals in this project are to develop extreme scaleusing the P-CG method.

AMR-CFD simulations and to establish performance

In two-phase flow simulations, a large density ratio

portability on the latest CPU and GPU platforms. For

between gas and liquid increases the contrast of each

this purpose, we need to use several supercomputing

matrix component. In addition, large-scale simulations

systems, which are based on state-of-the-art CPUs andhave multi-scale features such as thin gas-liquid

GPUs and have a large number of computing nodes.

interfaces and complex object boundaries. These

Such computing needs can be satisfied only by afeatures deteriorate the convergence property, and thus,

JHPCN joint research project.

the improvement of preconditioning is of criticalimportance.

4.

Outline of Research Achievements up to(3) Geometric Multigrid Preconditioner

FY2019The MG method is one of the most efficient

Not applicable.

preconditioners to reduce the computational cost andimprove the convergence property in multi-scale

5.

Details of FY2020 Research Achievements

problems. Since JUPITER is based on the block-5.1 JUPITER Codestructured Cartesian grid system, we can use the(1) Incompressible Two-Phase Flow Solver

geometric multigrid (GMG) method with a three-stage

We solve two phase flows using an incompressible fluid

V-cycle shown in Fig. 2. In the GMG method, anmodel based on the Navier-Stokes equation. The

approximate equation is solved at each MG level. The

Navier-Stokes equation consists of advection, diffusion,interpolation coefficients between MG levels are

and external force terms with an explicit time

calculated geometrically.

integration method. Since the divergence-free condition

The choice of the smoother at each MG level is veryon the velocity field is essential for mass conservation,

important in terms of the computational cost and thethe pressure Poisson equation is solved iteratively. As aconvergence property. The red and black SOR (RB-pressure-velocity correction algorithm, the simplified

SOR) method is an efficient smoother, which has highmarker and cell (SMAC) method is applied to JUPITER.

thread parallelism. However, it requires multiple

After solving the pressure Poisson equation, thevelocity field is modified with the updated pressure
gradient. The time evolutions of the above NavierStokes equation.
The dynamics of the gas and liquid phases aredescribed by a convection equation of the volume of
fluid (VOF) function. A physical property in a cell isdetermined by the linearly weighted average of theFig.2: The geometric multigrid preconditioner based on athree-stage V-cycle with fine, middle, and coarse MG
levels. The MG levels are defined within each leaf, inwhich the number of grids is reduced as 83, 43, and 23. A
fixed iteration number is used for smoothers at each level.

VOFs for the gas and liquid phases.
(2) Poisson Solver
The main computational cost comes from the pressure
Poisson equation discretized by the second order3

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2020

memory accesses for the red and black phases

Structured Cartesian Gridalternatively to eliminate the dependency between

JUPITER is written in C++ and CUDA. The connection

adjacent grids. In this research, we propose a cacheof the block-structured grid is managed by the forest-of-

reuse SOR (CR-SOR) method, which can effectively

octree data format, in which each array is based on the

reduce the residuals of local subdomain by multiple

Structure of Array (SoA) memory layout. The same

sub-iterations.

block-structure grid is used also in the MG method,(4) Block-Structured Cartesian Grid

where the offset index changes to 1/8 and 1/64 of the

The original JUPITER code was developed based on a

original resolution depending on the MG level.

Cartesian grid system, whose resolution was chosen to

It is important to assign CPU, GPU and A64FX

reproduce the shape of the VOF function or the gas-threads to grids in stencil computation. We implement aliquid interface with enough accuracy. However, there

kernel

was a lot of memory wastage for nuclear engineering“FOR_EACH1D_BLOCKIDX” corresponds to the

problems which contain complex non-fluid objects. In

leaves, and inner loop as “FOR_EACH3D” is assigned

the block-structured AMR grid system, grid resolution

to 83 grids in the leaf. Although the block-structured

is adaptively changed in each subdomain depending on

Cartesian grid system keeps continuous memory accessthe local scale of multi-scale solutions. Even without theas in the standard Cartesian grid system, special care haschange of grid resolution, the block-structured grid is

to be taken in the treatment of the leaves. Stencil

beneficial, when the boundary is not aligned to thecomputation for the block-structured grid requires the

Cartesian grid and/or complex structures exist inside theoffset indices of the neighbor leaves, which arecomputational domain. In this work, we use the block-normally assigned to the register memory. One leaf

structured grid with uniform grid resolution and reducerefers to the offset indices of 26 surrounding leaves, andthe grids outside the boundary and inside the object.

these indices can be shared by threads belonging to the

A data structure is very important from the

as

shown

in

Fig.

3.

Outer

loop

as

same leaf.

viewpoint of the computational efficiency on CPU/GPU
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:

platforms. A block-structured grid has an efficient data
structure, in which the whole domain is decomposed
into block units named “Leaf”. Here, the leaves areconnected by the forest-of-octree data structure. Since
the leaf consists of a regular Cartesian grid system with
n3 grids, continuous memory access is enabled. In thisresearch, the size of leaf is set to n=8 by taking accountof the above tradeoff.
The block-structured grid can be easily extended to
the MG method. It is possible to use the same forest-ofoctree structure as long as there is a grid in the leaf. We
adopt the MG method with a three-stage V-cycle, inwhich the number of grids at fine, middle, and coarse
MG levels are set to 83, 43, and 23, respectively.

#ifdef __CUDA_ARCH__#define FOR_EACH1D_BLOCKIDX(L, NL) \const auto L = blockIdx.x#define FOR_EACH3D(I,J,K, NX,NY,NZ) \const auto I = threadIdx.x \const auto J = threadIdx.y \const auto K = threadIdx.z#else#define FOR_EACH1D_BLOCKIDX(L, NL) \_Pragma(“omp parallel for”) \for(int L=0; L<NL; L++)#define FOR_EACH3D(I,J,K, NX,NY,NZ) \PRAGMA_FOR_SIMD \for(int K=0; K<NZ; K++) \PRAGMA_FOR_SIMD \for(int J=0; J<NY; J++) \PRAGMA_FOR_SIMD \for(int I=0; I<NX; I++)
#endifvoid cfd_function(...){
FOR_EACH1D_BLOCKIDX(l, num_amr_leaves){
FOR_EACH3D(i,j,k,NX_LEAF,NX_LEAF,NX_LEAF){// operations at stencil (i,j,k, l)}}}

Fig.3: Pseudocode for stencil computation to call CPU and
GPU instructions.

5.2 Implementation and Optimization(1) Implementation of JUPITER Code on Block4

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2020

Table 1: Hardware metrics of the TSUBAME and Flow
platforms.

(2) Cache Reuse SOR Method
The calculation speed of smoother is mainly limited bymemory read/write time rather than the cost of floatingpoint operations. The RB-SOR method requires global

Performance

memory access twice for the red and black phases,

TSUBAMETSUBAMEFLOW

NVIDIA

Broadwell

A64FX

P100

(2 sockets)(1 CMG)

4700

851.1

844.8

732

153.6

256

[GFlops]which causes performance degradation. To resolve thisBandwidth

issue, we propose a new CR-SOR method [Onodera etCompiler:

nvcc -O3

icpc -O3

mpiFCCpx

compiler-restrict-qopenmp-O3 -Kfast,

options-use_fast_-AVX2 -ipo

parallel, openmp,

math-restrictrestp=all[GByte/s]al., HPC Asia 2021]. The computational domain isdivided hierarchically in order to resolve dependencies
between leaves and grids inside the leaf. One is coarse

Table 2: Kernel performance analysis: Floating point
operation f [Flop/grid], Memory access b [Byte/grid],
Elapse time t [ns/grid], Roofline time tRL = f /F + b/B[ns/grid]. Here F is Peak performance [Flops], and B ismemory band-width [Byte/s]. (Note: The symbol * indicatesthe performance after optimization.)

domain decomposition for MPI parallelization (MPI
domain), and the other is fine block decomposition forcache optimization (cache block). In the outer iteration,
each MPI domain is processed for n times with a twocolor ordering. In the inner iteration, the RB-SOR
method within each cache block is computed for m

PCG kernels

times to accelerate the convergence.
f/b

5.3 Numerical Experiment(1) CG Kernel Performances on TSUBAME and

AXPY

FLOW
We evaluated the kernel performance of CG solvers in
hardware metrics of GPU and CPU on TSUBAME, and
A64FX on FLOW. Table 2 shows the kernel
performance analysis against the modified rooflinemodel [Shimokawabe et al., SC’10, 2010].

24

Inner

2/

product

16

SpMV

JUPITER on multiple platforms. Table 1 shows

3/

13 /
72

CR-SOR257 /(16) Lv.0

80

CR-SOR129 /(8) Lv.1

80

CR-SOR129 /(8) Lv.2

80

TSUBAMETSUBAMEFLOW

NVIDIA

Broadwell

A64FX

P100

(2 sockets)(1 CMG)

t

tRL/t

t

tRL/t

t

tRL/t

0.045

0.75

0.26

0.62

0.13

0.74

0.033

0.68

0.22

0.50

0.09

0.72

0.141

0.72

1.19

0.41

7.45

0.04

*0.51

*0.58

0.447

0.37

6.03

0.14

72.7

0.01

0.368

0.37

3.75

0.18

26.4

0.02

0.472

0.29

3.09

0.22

25.3

0.02

AXPY and inner product kernels are memory
intensive kernels without neighbor stencil access. Since

in a loop, it was necessary to perform manual inline

the block-structured data format enables continuous

expansion for those functions and to replace thememory access in the leaf, these kernels achieved ideal

conditional branch with integer flags. Thanks to thesesustained performance with the performance ratio

optimization techniques, finally, the performance on

against the roofline model tR/t > 0.5.

A64FX became comparable that on CPU and GPU.

SpMV kernel is a very important kernel, because

CR-SOR kernel is computationally intensive kernel

SOR kernels also have similar memory access. The

on each MG level. Memory accesses are almost theperformances on GPU and CPU are reasonable, butsame as SpMV kernel, but this kernel computes

there is a significant performance degradation on

multiple SpMV processes on a local cache. The

A64FX. According to the optimization report, theperformances on GPU and CPU were reasonable, tR/t >innermost loop was not SIMD-optimized due to inline

0.35 on GPU and tR/t > 0.15 on CPU. On the other hand,functions and a conditional branch in the loop. Since

the performance on A64FX was quite low, tR/t < 0.02,

Fujitsu C++ compiler cannot optimize any function callwhile we expected the similar sustained performance as
5

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2020

SpMV kernel. Although we investigated this issue withour collaborators from Fujitsu, the performancedegradation was not resolved even with the above
manual optimization. Possible causes for this issue may
be deeper nested loops due to inner iterations and the(A)two-color ordering or very tight scheduling for on-cachecomputing.

(B)

Fig.4: (A) shows the computational domain for a 8×8 fuelbundle. The size of the whole domain is 104 mm×104 mm
×312 mm. (B) shows the block-structured grid. Blue andred leaves show fluid and boundary regions, which areallocated in the simulation. Leaves inside the pins (whitearea) are not allocated.

(2) Convergence Histories of MG-CG Method
The numerical experiments were conducted for twophase flows in a fuel bundle of a nuclear reactor shownin Fig. 4. Since computational grids inside fuel pins are
removed, the total number of grids is reduced to
2.26×109, which is about 70% of the original uniform
Cartesian grid system.
Figure 5 shows convergence histories of the P-CG
solver and the MG-CG solvers with the RB-SOR
smoother (MG-CG) and with the CR-SOR smoother(CRMG-CG). Although the P-CG solver required morethan 1,300 iterations until convergence, the number of
iterations in the MG-CG solver was 194, which was lessFig.5: Convergence histories of the P-CG, MG-CG, and
CRMG-CG solvers. A large matrix from a two phase flow
problem in a 8×8 fuel bundle is computed using 64 GPUs.
The P-CG, MG-CG, and CRMG-CG solvers arerespectively converged with 1,306, 194, and 114 iterations.
The calculation time until convergence of the P-CG, MGCG, and CRMG-CG solvers are respectively 446, 142, and76 seconds.

than 15% of the P-CG solver. The CRMG-CG solver
showed faster convergence with 114 iterations, which
was less than 9% of the P-CG solver.
(3) Strong Scaling Test on GPU
The strong scaling of JUPITER with the MG-CG and
CRMG-CG solvers is summarized in Fig. 6. In thestrong scaling test, we use 64, 128, 192, and 256 GPUs.
The labels of “Precondition”, “CG”, and “NS” indicate
the preconditioner, the CG method including the innerproduct, and the Navier-Stokes solver with the explicit
time integration, respectively.
The breakdown of computational time on 64 GPUs
shows that Precondition accounts for more than 90% of
the total time, the remaining cost comes from CG, andNS is negligible. A similar breakdown was observed up
to 256 GPUs. In the MG-CG and CRMG-CG solvers,Fig.6: Strong scaling tests of JUPITER with the MG-CG
and CRMG-CG solvers up to 256 GPUs. The
computational condition is the same as Fig. 9. It is notedthat the total time of JUPITER with the P-CG solver is447.5 second on 64 GPUs.

global collective communication in the CG method isstill negligible, and halo data communication gives adominant communication cost. The performance gainsof the CRMG-CG solver from 64 GPUs are 1.5×, 1.9
6

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2020

×, and 2.1× on 128, 192, and 256 GPUs, respectively.

Conjugate Gradient Solver on Block-Structured

It is noted that the calculation speed of the CRMG-CG

Cartesian Grid”, Proceedings of HPC Asia 2021,

solver on 256 GPUs is 12 times faster than that of the2021/01, 査読有り

P-CG solver on 64 GPUs.



N. Onodera, Y. Idomura, Y. Ali, S. Yamashita,
T. Shimokawabe, T. Aoki, “GPU-acceleration of

6.

Progress during FY2020 and Future Prospects

locally mesh allocated two phase flow solver forThis project developed a new MG-CG solver tonuclearaccelerate a three-dimensional two-phase flow CFD

SNA+MC2020, 2020/05, 査読有り

reactors”,

Proceedings

code JUPITER. The MG-CG solver is designed to(3)International conference Papers (Non-refereed)achieve high performance on the block-structured

N.

Cartesian grid.

Onodera,Y.

Idomura,Y.

Asahi,ofY.

Hasegawa, S. Yamashita, T. Shimokawabe, T.

We measured the kernel performance of the MG-CG

Aoki, “Multigrid Poisson Solver for a Block-solver on CPU and GPU on TSUBAME, and A64FX on

Structured Adaptive Mesh Refinement Method

FLOW. Although AXPY, inner product, and SpMVon CPU and GPU supercomputers”, Compsafe

kernels achieved reasonable performances on all

2020, 2020/12, 査読無し

platforms, CR-SOR kernels showed a significant(4)

performance degradations on A64FX. In the futurerefereed)

project, this issue will be investigated on FLOW and
BDEC

under

collaboration

withPresentations at domestic conference (Non小野寺 直幸, 井戸村 泰宏, 長谷川雄太, 山下晋,

Information

下川辺 隆史, 青木 尊之, “ブロック型適合細分化

Technology Center, Nagoya Univ. and Information

格子を用いた気液二相流体解析手法の開発”, 原Technology Center, Univ. Tokyo.

子力学会 2021 年春の年会, 2021/3, 査読無し

We discussed the convergence property for two-

小野寺 直幸, 井戸村 泰宏, 朝比 祐一, 下川辺隆史, 青木 尊之, “ブロック型適合細分化格子での Poisson 解法の GPU・CPU・ARM プロセッサに対する性能測定”, 第 34 回数値流体力学シンポ
ジウム, 2020/12, 査読無し



小野寺 直幸, 井戸村 泰宏, ユスフ アリ, 下川辺隆史, 青木 尊之, “ブロック型適合細分化格子での Poisson 解法の GPU 高速化”, 計算工学講演会
論文集, 2020/06, 査読なし

phase flows in a fuel bundle. The number of iterationsof the MG-CG and CRMG-CG methods wererespectively less than 15% and 9% compared to the PCG method.
The strong scaling of JUPITER was measured from64 to 256 GPUs. The performance gain of the CRMGCG solver from 64 GPUs to 256 GPUs was 2.1×. We
conclude that the redesigned JUPITER code is highly
efficient and enables large-scale simulations of twophase flows on GPU based supercomputers.

7.

List of Publications and Presentations(1)

Journal Papers (Refereed)(2)

ProceedingsofInternational

Conferences(Refereed)

N. Onodera, Y. Idomura, Y. Hasegawa, S.
Yamashita, T. Shimokawabe, T. Aoki, “GPU
AccelerationofMultigrid

Preconditioned

7