2020年度第1回計算科学フォーラム
2020年11月24日（火）15：00～17：00

高性能計算物理シミュレーションと人工知能の融合の試み-地震シミュレーションを例に市村強東京大学地震研究所計算地球科学研究センター
東京大学工学系研究科社会基盤学専攻
概要：地震問題の中には，対象が巨大かつ高分解能が求められるため大規模物理シミュレーションが必要となるものがある．本講演では，このような大規模物理シミュレーションを可能とするために開発を続けている高性能計算によるマトリクスソルバーを用いた非構造低次有限要素解析とそのシミュレーション能力をさらに高めることを目指した物理シミュレーションと人工知能の融合の試みについて，いくつかのスーパーコンピュータ上で実施した大規模地震シミュレーションとともに紹介する．

outline• ポ ス ト「 京 」重 点 課 題 3 「 地 震 ・ 津 波 に よ る 複 合 災 害 の 統 合
的 予 測 シ ス テ ム の 構 築 」で進められ，また，「『富岳』成果創出加速プログラム：大規模数値シミュレーションによる地震発生から地震動・地盤増幅評価までの統合的予測システムの構築とその社会実装」で現在進められている以下について概説する• High Performance Computing (HPC) による物理シミュレーションの高度化• HPC物理シミュレーションの高度化と人工知能の融合へ向けて• HPC物理シミュレーションによる人工知能の高度化• 人工知能によるHPC物理シミュレーションの高度化• 人工知能アーキテクチャー？によるHPC物理シミュレーションの高度化• 微分方程式を学習した人工知能によるHPC物理シミュレーションの高度化

High Performance Computingによる物理シミュレーションの高度化
京コンピュータを用いた首都圏の大規模地震シミュレーション
Tsuyoshi Ichimura, Kohei Fujita, Seizo Tanaka, Muneo Hori, Maddegedara Lalith, Yoshihisa Shizawa, and Hiroshi
Kobayashi, Physics-based urban earthquake simulation enhanced by 10.7 BlnDOF x 30 K time-step unstructured
FE non-linear seismic wave simulation, SC14: International Conference for High Performance Computing,Networking, Storage and Analysis, 15-26, 2014.
Tsuyoshi Ichimura, Kohei Fujita, Pher Errol Balde Quinay, Lalith Maddegedara, Muneo Hori, Seizo Tanaka,Yoshihisa Shizawa, Hiroshi Kobayashi and Kazuo Minami, Implicit Nonlinear Wave Simulation with 1.08T DOF
and 0.270T Unstructured Finite Elements to Enhance Comprehensive Earthquake Simulation, SC15: International
Conference for High Performance Computing, Networking, Storage and Analysis, Article No. 4, 2015.

解析対象6

背景１

地表面の揺れSoft soil layers

入力基盤地震動<100m• 地震対策用シナリオ（都市・構造物の被害想定等）構築のためには，地震時の地表面の揺れ分布の情報が面的に必要• 地震の揺れは場所によって大きく異なる．

• 地震時に地盤材料は強い非線形性を示す．
• 特に，ローカルな浅い３D地盤構造が効く（深さ100ｍくらいまでの軟らかい層）

局所的な地盤構造が地震の揺れ分布に効く一例

軟らかい層

＜100ｍ

1920m
1696m

地表面での地面の揺れの強さ分布
・＜50ｍ程度の分解能で揺れが大きく異なる背景２• 10ｍオーダーの起伏を有する３D地盤構造のデータ整備がすすめられている• 例えば，左図四角内の地盤構造

地表面標高

5 miles (8.04km)

1層目と2層目の層境標高

2層目と3層目の層境標高
全国電子地盤図より作成

地盤増幅解析のための非線形有限要素解析• 低次非構造要素を用いた陰的時間積分非線形波動有限要素解析が必須• 地盤構造の幾何形状を正確にモデル化・応力フリーの境界条件を解析的に満足• 高い不均質性・強い非線形性を解析• 安定性の観点から陰的時間積分が必要
Matrix (components
Unknown vector
changes every time step) (1000億自由度)

Known vector
この離散化された非線形波動方程式を6000時間ステップ解く
陽的時間積分に比べ格段に難しい∵方程式を解く必要＆並列計算のglobal communication

大規模有限要素解析手法
Solving preconditioning matrix
CG loop

Inner coarse loop

Solvingpreconditioning
matrixSecond orderedtetrahedron

Computations ofouter loop

Equation to be solved(double precision)

Outer loop

Solve system roughly using CG solver(single precision)
Use as initial solution
Inner fine loop

Lineartetrahedron

Solve system roughly using CG solver(single precision)

Use for preconditioner of outer loop

• Solve preconditioning matrix roughly to reduce number of CG loops

Second orderedtetrahedron• Use multi-grid method to reduce cost of preconditioner• Use single precision in preconditioner to reduce computation & communication• Structuring unstructured mesh

• Omit costly geometry computation part of unstructured elements

大規模有限要素解析手法の効率• Size-up scalability: 96.6% efficiency from 9,216 cores to full K computer(663,552 cores)• Enables 18.6% of peak (=1.97 PFLOPS) & 1.08 trillion DOF @ full K computer• 8.6 times faster than conventional solver (CG + Element-by-Element + simple
preconditioner (block diagonal Jacobi preconditioner)• Speed-up scalability: 76% efficiency for 9,216 ➔ 294,912 cores• Very good scaling considering non-uniform mesh is partitioned using METIS• Similar scalability can be attained for practical problem

Normalized
elapsed
time

29216
1
1/2
1/4
1/8
1/16
1/32

294912 663552 # of coresSize-up scalability

断層-都市-社会の一気通貫地震シミュレーション• 首都圏の山手線内を対象とした仮想首都直下地震のシミュレーション（京コンピュータのmax能力で）

5 miles (8.04km)

断層から地表までの地震動シミュレーション• 10Hzまで精度保証した非構造有限要素による地震動シミュレーション• 18 billion の低次非構造有限要素＆56 billion 自由度による陰的時間積分

UNECity simulation domain
Fault

Velocity Norm (cm/s)

10.0
5.0
T=15 s

T=16 s

T=17 s0.0

16City-to-social simulation

Visualized by Leonel Aguilar (University of Tokyo)

First ever practical fault-to-city-to-social simulationaccomplished• All phases of earthquake fault-to-city-to-social simulation
simulated at state-of-art level• Enables huge leap in improving reliability of earthquake
response estimation• HPC can play an essential role for saving lives and society

Fault-to-city simulation

City simulation

Social simulation 18

Equation-based modelingとdata analytics
Equation-based modeling

解微分方程式等
数式モデル

現象

差分法・有限要素法等
数値シミュレーション

数値解

データ

深層学習等
data analytics

予測モデル

Equation-based modelingとdata analyticsの融合によりより多くの問題が解けるようになる?

HPC物理シミュレーションの高度化と人工知能の融合へ向けて
HPC物理シミュレーションによる人工知能の高度化
Tsuyoshi Ichimura, Kohei Fujita, Takuma Yamaguchi, Muneo Hori,
Maddegedara Lalith, Naonori Ueda, AI with Super-Computed Data for Monte
Carlo Earthquake Hazard Classification, SC17: International Conference forHigh Performance Computing, Networking, Storage and Analysis, 2017.

地盤構造の曖昧さ地表面の揺れ
Soft soil layers

入力基盤地震動?

<100m• 地表面は直接計測可能だが，地盤内部は直接計測はできな
いため，精度があまり高くな
く，あいまいさあり．
• あいまいさを考慮するための対処が必要• 強い非線形性があるので，一番良いのはモンテカルロシミュ
レーション（先ほどの解析を1000回実施するか？）

超大規模解析による超大規模データ

10250m
3 soil layers

9250m• 京コンピュータ全系(663552 CPUコア)を使用上限（8時間）用いて実施• 最小要素サイズ（非構造四面体二次要素）: 1m• 1333億自由度で6600時間ステップ• 情報量: 1333億自由度 X 倍精度 X
6600 時間ステップ = 6600 TB

地震の揺れマップ 物理シミュレーション vs AI

By Supercomputing

By AI with Super-computed data

両者は良く一致している∵人工知能の正答率が高くなるように学習させる方法論を開発提案したため

人工知能+Monte-Carloによる地震の揺れマップ

mean

3 X Standard deviation

1000 trialsによる曖昧さを考慮した地震の揺れマップ

・地震の揺れマップに地盤構造のあいまいさが及ぼす影響は場所によって大きく異なる・meanだけでなく，standard deviationも場所によって大きく異なることがわかる

HPC物理シミュレーションの高度化と人工知能の融合へ向けて
人工知能によるHPC物理シミュレーションの高度化
Tsuyoshi Ichimura, Kohei Fujita, Takuma Yamaguchi, Akira Naruse, Jack C. Wells,
Thomas C. Schulthess, Tjerk P. Straatsma, Christopher J. Zimmer, Maxime Martinasso,
Kengo Nakajima, Muneo Hori, Lalith Maddegedara, A Fast Scalable Implicit Solver for
Nonlinear Time-Evolution Earthquake City Problem on Low-Ordered Unstructured
Finite Elements with Artificial Intelligence and Transprecision Computing, SC18:International Conference for High Performance Computing, Networking, Storage andAnalysis, Article No. 49, 2018.

Equation based modelingとdata analytics• Equation based modeling• 高精度, 解析コスト高• Data analytics• 高速に評価可能， Equation based modelingと比較して精度低Phenomena

Data analytics

Equation based modelingEquation based modelingとdata analyticsの融合へ• 第一段階: equation based modelingによって生成されたデータをdata analytics に用いる
• 大規模なequation based modelingにより膨大な高品質のデータを生成可能• 実験が難しい/不可能な現象の疑似観測データを生成• 地震の揺れの強度を推定するAIをこの手法で開発した
SC17
Phenomena

Data analytics(with betterprediction)

Simulated data Equation basedmodelingfor training

•
•
•
•

SC14: equation based modeling
SC15: equation based modeling
SC16: equation based modeling
SC17: equation based modeling for AI

Equation based modelingとdata analyticsの融合へ• 第二段階：さらなる融合へ．AIにより equation basedmodelingを高速化・解析コスト低減へ
SC18

Dataanalytics

Phenomena

AI for accelerating
equation basedsolver

Equation basedmodeling(25-fold speedupfrom without AI)

••
•
•
•

SC14: equation based modeling
SC15: equation based modeling
SC16: equation based modeling
SC17: equation based modeling for AI
SC18: AI for equation based modeling人工知能により系を学習し，方程式求解を高速化• 微分方程式の持っている情報を用いる• 方程式の解きにくい部分の情報を支配方程式がもっている• 解きにくい部分をAIで抽出して，優先して解いてしまう

Phenomena

Dataanalytics

Equation based modeling
支配方程式（微分方程式）
離散化

Ax=fSummit上での性能（Weak scaling）
500
400

349.8

374.6

380.2

415.1

454.2

311.7

327.3

75.8

77.6

80.4

82.9

84.3

83.7

90

100.4

# of GPUs

288
288

576
576

1152
1152

2304
2304

4608
4608

6144
6144

12288
12288

24576
24576

# of node

48

96

192

384

768

1024

2048

4196

DOF

3.5×109

7×109

14×10928×10956×10975×109

151×109

302×109

SC18’s efficiency toFP64 peak

19.5%

SC18MOTHRASC14GAMERA

Elapsed time (s)

302.5

300
200
100014.7%

人工知能によるHPC物理シミュレーション on Summita) Overview of citymodelb) Location of underground structurec) Close up view of city modeld) Displacement response of citye) Displacement responseof underground structureHPC物理シミュレーションの高度化と人工知能の融合へ向けて
人工知能アーキテクチャー？によるHPC物理シミュレーションの高度化
Tsuyoshi Ichimura, Kohei Fujita, Takuma Yamaguchi, Akira Naruse, Jack C. Wells, C. J.
Zimmer, T. P. Straatsma, Takane Hori, Simone Puel, Thorsen W. Becker, Muneo Hori, andNaonori Ueda, 416-PFLOPS Fast Scalable Implicit Solver on Low-Ordered Unstructured
Finite Elements Accelerated by 1.10-ExaFLOPS Kernel with Reformulated AI-LikeAlgorithm: For Equation-Based Earthquake Modeling, SC19: International Conference forHigh Performance Computing, Networking, Storage and Analysis, 2019.

人工知能用？アーキテクチャの活用• 人工知能用？に開発された低精度・高速アクセラレータをHPC物理シミュレーションで活用できないか• 例えば，Tensor Core on NVIDIA Volta GPUは，FP64演算に比べて16倍高いpeak FLOPS：FP64 7.8 TFLOPS, Tensor Core 125 TFLOPS• 一方で，低精度演算，局所性の高い演算である必要• 定型の行列積などの制約も．．．

• 人工知能用？アーキテクチャを活用できるように，従来のHPC物理シミュレーションとは異なる上記の制約にあうアルゴリズムを開発• 一見すると非効率？

Tensor Coreを用いた物理シミュレーション

前処理部にTensor Coreにあわせて構築した構造格子的演算を導入
・FP16の16*16の行列積ベースの演算に再構築
・GPU内のデータ移動コストに配慮，低精度通信導入など計算機構を踏まえたアルゴリズム開発によりtime to solutionを改善．Flop countではなく，各計算機構での実効Flop countの観点が重要• 地殻変動解析用の地殻モデルを模した
1.67×1012自由度の問題の二層モデルを，Summit 4544 計算ノード (27,264 NVIDIAV100 GPUs) で解析• 提案手法では演算量が増加（主に構造格子的演算部）．しかし，Tensor Core が行列積を高速化することで，SC14Solverに比べて3.89倍の高速化を実現

人工知能アーキテクチャー？によるHPC物理シミュレーション on Summit• Cascadia沈み込み帯における断層すべりによる地殻変動の解析• 1944 km×2646 km×480 km 領域の1.49×1010 DOF の地殻モデル• Summitの162計算ノードを用いて23sで計算可能• state-of-the-art solver on Piz Daintと比べて75倍の高速化を実現• 416-PFLOPS Fast Scalable Implicit Solver on Low-Ordered Unstructured Finite Elements
Accelerated by 1.10-ExaFLOPS Kernel

Crust/mantle
interfacegeometry
Pacific
plate
Maximum principal
stress distribution

NorthAmerican plate

HPC物理シミュレーションの高度化と人工知能の融合へ向けて微分方程式を学習した人工知能によるHPC物理シミュレーションの高度化
Tsuyoshi Ichimura, Kohei Fujita, Muneo Hori, Lalith Maddegedara, Naonori
Ueda, and Yuma Kikuchi, A Fast Scalable Iterative Implicit Solver with
Green’s function-based Neural Networks, 2020 IEEE/ACM 11th Workshop
on Latest Advances in Scalable Algorithms for Large-Scale Systems(ScalA)@SC20: International Conference for High Performance Computing,Networking, Storage and Analysis, 2020.

概要• Based on the Green’s functions that reflect mathematical properties of partial
differential equations (PDE), we developed a novel preconditioner usingneural networks (NNs) with high accuracy and small computational cost for
improving the convergence property of an iterative implicit solver.
• As the dense and uniform computation involved in NNs are more efficientthan that of the conventional PDE solver schemes, we could solve the time
evolution of a highly heterogeneous problem in 4.26-fold shorter timecompared to a typical PDE solver. The method is also suitable for use with
low-precision arithmetic in NNs as the accuracy of the final solution isguaranteed. The localized property of NNs enable high scalability for solving
large problems: a 5.57 trillion degrees-of-freedom problem was solved with
98.1% weak scalability up to 55,296 Fugaku nodes with 30.3 PFLOPS (16.2%of FP64 peak performance). This method can be used in various PDE-based
simulations and has potential to make broad ripple effects in various fields.

Green’s function-based Neural Networks
を用いたimplicit matrix solver
非線形時刻歴発展問題の以下の求解を考える．

g(x, s)をNNを用いて構築することとする

ここでAとfは時刻n-step目に依存して変化する．

非線形波動場問題のGreen functionとGreen’s function-based Neural Networks

一方で，PDEにおいて，が期待される．ここで，(L(x): linear differentialoperator, b(x): known distribution, a(x): responseof this system, x: an arbitrary point).
g(x, s)がGreen’s functionであり，これを前処理として用いることを考える．
一見良さそうだが， g(x, s)の計算コストが大きすぎ使えない．．．

一般的なCGの前処理に用いた例

Distribution of the Young’s modulus of the target head model

Convergence speeds of solvers for the head model measured on one
Oakbridge-CX node. The convergence ratio (ratio of change in relative
error per iteration) was greatly improved by the use of the GF-based
NNs in each refinement iteration in the proposed method. The simple
scheme using GF-based NNs for initial solution of PDE solver scheme(CGBJ) could not utilize this improvement in the convergence ratio.

※Numbers in brackets indicate ratio to FP64 peak.

GF-based NNはsolverの収束率を改善可能だが，実は演算量自体は増加．
しかし，連続アクセス演算卓越なので高速に処理可能で，また，GF-based NNはポイ
ントワイズな前処理なのでscaleしやすい

Weak scaling results of proposed solver on Fugaku: a 5.57 trillion degreesof-freedom problem was solved with 98.1% weak scalability from 864 up to55,296 Fugaku nodes with 30.3 PFLOPS (16.2% of FP64 peak performance).

まとめ• ポ ス ト「 京 」重 点 課 題 3 「 地 震 ・ 津 波 に よ る 複 合 災 害 の 統 合 的 予 測シ ス テ ム の 構 築 」で進められ，また，「『富岳』成果創出加速プログラム：大規模数値シミュレーションによる地震発生から地震動・地盤増幅評価までの統合的予測システムの構築とその社会実装」で現在進められている以下について概説した• High Performance Computing (HPC) による物理シミュレーションの高度化• HPC物理シミュレーションの高度化と人工知能の融合へ向けて• 物理シミュレーションと人工知能の融合は他にも様々な方向性が考えられる．さらなる評価手法全体の性能向上が期待される• 従来の物理シミュレーションのアルゴリズムを人工知能的に改良• Flop countだけではなく，計算機構に沿った実効Flop countの観点からより強力なアルゴリズムを設計が可能に など• Society 5.0へ向け，IoTや５Gを踏まえた観測と解析の融合（Big Data & ExtremeComputingなど）への展開

以上のような検討により高効率な大規模都市解析が可能になりつつある（人工知能によるHPC物理シ
ミュレーション@富岳の例）

Computation example of 11,002,859,706 degrees of freedom; 2,593,545,132 element model of fully
coupled ground and aboveground/underground structure earthquake shaking analysis using IRIS on
98,304 CPU cores (2,048 nodes) of Fugaku: The 2 x 2 km domain was modeled with 0.5-m tetrahedral
elements with three soil layers; 9323 buildings; and an underground complex with 645 tenants. 44

Acknowledgement• We acknowledge supports from Program for Promoting Researcheson the Supercomputer Fugaku (Large-scale numerical simulationof earthquake generation, wave propagation and soil amplification:
hp200126) and Post K computer project (Priority Issue 3 Development of integrated simulation systems for hazards anddisasters induced by earthquakes and tsunamis: proposal numbers
hp180217, hp190177), and Japan Society for the Promotion ofScience (18H05239, 26249066, 25220908, and 17K14719). Our
results were obtained using Fugaku and K-computer at Center for
Computational Science, RIKEN and Oakbridge-CX at Information
Technology Center, The University of Tokyo and the Summit at OakRidge Leadership Computing Facility, a US Department of Energy,
Office of Science User Facility at Oak Ridge National Laboratory.