2022/10/27

データ・コンピューティング
による高速解析- 地震シミュレーションを例に
市村強（東京大学・地震研究所）

解けない問題が多々ある
解析手法開発解析能力をさらに向上させるニーズ（従来法の枠からはみ出した？新たな方法へ．．．）

新たな「先端的解析手法」の開発研究課題の提案へ

先端的解析手法の開発
科研費研究：基盤S，萌芽など相補的

新たな「超大規模問題への展開」の開発研究課題の提案へ

先端的解析手法の超大規模問題への展開「富岳」成果創出加速プログラム：「大規模数値シミュレーションによる地震発生か
ら地震動・地盤増幅評価までの統合的予測システムの構築とその社会実装」※富岳課題概要説明は次頁から

実現象の解明・災害軽減へ向けた各PJ

2022/10/27

本課題の位置づけ・最終目標

by courtesy of
Dr. Takane Hori

• 「大規模数値シミュレーションによる地震発生から地震動・地盤増幅評価までの統合的予測システムの構築とその社会実装」• 応募領域:②国民の生命・財産を守る取組の強化

国の被害想定のための地震・ 津波災害予測における ア プ リ ケーシ ョ ン 活用・ 重点課題から 開発し て き た大規模有限要素を 用いた計算コ ード が、 国の被害想定のも と になる 長周期地震動計算や津波初期水位計算で 用いら れる よ う ア プ リ ケーシ ョ ン の改良を 国と 連携し て 実施する 。
・ 建築・ 土木系企業等が、 国の想定と 同等の計算がで き る 仕組みを 構築する 。

地震に関する 災害被害予測のための大規模ア プ リ ケーシ ョ ン の開発
・ 「 富岳」 の性能を 引き 出すよ う に、 計算科学・ 計算機科学の最先端技術を 駆使し て 、 地震に関する 災害被害予測のための大規模ア プ リ ケーシ ョ ン を 改良する 。

地震災害の定量的評価を超大規模計算で実施するための統合的予測システムを構築して、国の被害想定における地震動・地盤増幅評価の過程に実装し、政策判断や実務等に活用される基礎を築く。

実施体制・関係機関

by courtesy of
Dr. Takane Hori

共通基盤：超大規模高速FEMソルバー
HPC・AI技術を駆使し、「富岳」の性能を引き出すFEMソルバ―の開発

カーネルコード開発：東大地震研（市村・藤田）

テーマ１：長周期地震動

テーマ２：地殻変動

テーマ３：浅層地盤増幅

地震波伝播、長周期地震動手法の開発と適用

地震発生、地殻変動手法の開発と適用

浅層地盤における非線形地盤
増幅手法の開発と適用企業における
長周期地震動計算：
大崎総研・
小堀鐸二研・
大成建設

国の長周期地
震動計算への適用：
内閣府防災担当・東大地震研
災害科学

各テーマ向けコード改変
大規模計算
実計算用機能・実行環境整備

企業との共同研究

企業における非線形地盤増幅
計算：東電設計・東京ガス

プロジェクト
の統括

国の被害想定への適用：
内閣府防災担当

東大地震研（市村・藤田）海洋研究開発機構（堀・縣）

国・企業・サイエンスコミュニティとの連携

計算結果のサイエンスコミュニティでの活用
東大理固体地球・科研費新学術「スロー地震」 ・地震調査研
究推進本部・防災科研・文科省委託事業「防災南海P」

国の津波被害想定への適用：
内閣府防災担当

計算結果の都市被害予測への活用
理研R-CCS総合防災・防災科研代表機関 協力機関 連携機関 関係事業・機関○代表機関:

海洋研究開発機構（プロジェクト代表:堀高峰）○協力機関:

東京大学地震研究所（プロジェクト副代表:市村強）○連携機関:

内閣府（防災担当）、東京大学大学院理学系研究科、理研R-CCS
株式会社 大崎総合研究所、株式会社小堀鐸二研究所、大成建設株式会社、東京ガス株式会社、東電設計株式会社

2022/10/27

これまでの成果by courtesy of
Dr. Takane Hori

国の被害想定のための地震・津波災害予測でのアプリケーション活用• 「相模トラフ沿いの巨大地震に伴う⾧周期地震動による被害予測の高度化」
内閣府政策統括官（防災担当）の政策枠での実施• 工学的基盤での周期2秒までの線形波動場計算を関東周辺域で実施

地震に関する災害被害予測のための大規模アプリケーションの開発• 関東平野を含む広域モデルで、地下深部の地震発生から構造物の揺れまでのシームレスな計算を「富岳」全系で実施• 関東平野周辺から地下＆地上構造物まで一気通貫の地震動計算を実現開発した計算アプリケーション群を実務で活用するための環境整備• 国の⾧周期地震動計算と同等のアプリのユーザーインターフェースを整備・改良し、現実的な地下構造を対象とした問題に連携企業の方々が適用

背景• データの蓄積により，多くの良好な三次元モデルが提案されてきている（地殻構造モデル，地盤モデ
ル等の理～工スケールで）• 対象の挙動は，静的/動的，非線形/線形の応答としてモデル化される

理～工の対象とする問題全体像

数理問題の観点からみると，見通しがつきやすい（レングスケー
ルや分解能が違う位で，理と工をそこまで区別する必要はない）

・これらの順・逆解析をする問題に帰着
・問題によって「構成則」や「感覚」は違うが．

2022/10/27

目的（応力フリーの）境界条件や幾何形状が解に強く影響を及ぼす数理問題なので，有限要素法を使いたい．これが出来れば，理も工も共通の解析基盤にのる？
Continental
plate
Pacific
plate
Philippine
sea plate

日本列島規模784km

y

町内会？規模1184kmzx

観測データ＋解析→virtual 地球 on computer

400km

東日本規模

地下高速道路トンネル

理も工も関係なく，観測と解析をスムーズにつなげる（自分以外の分野も分かりやすくなる）

Society 5.0 のためのBig data & extreme computing
観測データ（都市観測データなど）
計測

5GやIoTにけん引されるDXにより従来より巨大化→Supercomputingがあればさらなる深化が期待される解析結果を受けた検討・判断など
実空間

都市のデジタルツイン

分析・推定・予測

観測データに整合する計算機上に構築された双子（ex) デジタルツイン東京）：
データ群とシミュレーション群の塊

データフロー
解析
ストレージ

デジタル空間Big Data & Extreme Computing System
フローからの必要情報及びパラメータ等の解析結果をアーカイビング

多数の観測機器から得られたデータを蓄積し，大規模・多種多様なシミュレーションを駆使して過去・現在・未来を評価可能な都市・社会・地球の次世代デジタルツインによる未来予測・設計制御へ

2022/10/27

Seismic structural response analysis
with high fidelity modela) overview of the plant model with surrounding
two-layered ground modeled with 11,321,249,889
elements and 49,064,764,344 degrees of freedom;b) close-up; c) complex geometry modeled with
16-mm sized second-order tetrahedral elements;d) displacement response to ground shaking; ande) stress concentration, which can be evaluated forall parts of the model at high-resolution forengineering purposes.
Tsuyoshi Ichimura, Kohei Fujita, Masashi Horikoshi, Larry
Meadows, Kengo Nakajima, Takuma Yamaguchi, Kentaro
Koyama, Hikaru Inoue, Akira Naruse, Keisuke Katsushima,
Muneo Hori, Maddegedara Lalith, A Fast Scalable Implicit
Solver with Concentrated Computation for Nonlinear Timeevolution Problems on Low-order Unstructured Finite
Elements, IEEE International Parallel and Distributed
Processing Symposium, 2018.
9

大規模地震シミュレーション
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.
10

2022/10/27

地震解析のための非線形有限要素解析• 低次非構造要素を用いた陰的時間積分非線形波動有限要素解析が必須• 構造の幾何形状を正確にモデル化・応力フリーの境界条件を解析的に満足• 高い不均質性・強い非線形性を解析「解の収束性確認」の重要さ• 安定性の観点から陰的時間積分が必要
Matrix (componentschanges every time step)

Unknown vector(1000億自由度)

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

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

Inner coarse loop

Solvingpreconditioning
matrixSecond ordered
tetrahedron

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 ordered
tetrahedron• 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

2022/10/27

超並列大規模有限要素解析手法の効率-富岳，京，Summit，OFP等で• 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

2
1

9216

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

294912 663552# of coresSize-up scalability

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

5 miles (8.04km)

2022/10/27

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

N
E

City simulation domain
Fault

T=15 s

T=16 sVelocity Norm(cm/s)

T=17 s16

2022/10/27

City-to-social simulation

Visualized by Leonel Aguilar (University of Tokyo)

データ・コンピューティング：
系の性質を反映したデータによるコン
ピューティングの加速微分方程式を学習した人工知能による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 andAnalysis, 2020.
Kohei Fujita, Yuma Kikuchi, Tsuyoshi Ichimura, Muneo Hori, Lalith Maddegedara, Naonori Ueda, GPU
porting of scalable implicit solver with Green’s function-based neural networks by OpenACC,
WACCPD 2021: Eighth Workshop on Accelerator Programming using Directives@SC21.

2022/10/27

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

g(x, s)をNNを用いて構築することとする
非線形波動場問題のGreen functionと

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

Green’s function-based Neural Networks

において，が期待される．ここで，(L(x): linear differential operator,b(x): known distribution, a(x): response of 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.

2022/10/27

時刻歴解析結果（uの分布）
上段：従来のアプローチ（ブロックヤコビ＋CG）
中段：人工知能のみによる
求解（NNベースG）
下段：提案手法（人工知能前処理＋CG）
・高精度な人工知能を構築出来ているが，これだけでは精度が劣化していく
・CGと組み合わせることで，精度劣化なく解析可・人工知能解析の精度保証をCGにより実現
・人工知能によるCGの高速化の二つの側面※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 to
55,296 Fugaku nodes with 30.3 PFLOPS (16.2% of FP64 peak performance).

2022/10/27

GPU(A100) porting with OpenACC
ABCI A100 1 computer node
Total Hardware peak

CPUIntel Xeon Platinum
8360Y (2.40 GHz, 36
Cores) × 2

5.529 TFLOPSGPUNVIDIA A100 NVLink×8

77.6 TFLOPS

Time-to-solution of the whole solver (proposed method)
CGBJ@CPU
Proposed@CPU
Proposed@GPU(Fix the problem size to 256 × 256 × 512 elements per GPU)

データ・コンピューティング：
系の性質を反映したデータの逐次学習によるコンピューティングの加速
Tsuyoshi Ichimura, Kohei Fujita, Kentaro Koyama, Ryota Kusakabe, Yuma Kikuchi, Takane Hori, Muneo
Hori, Lalith Maddegedara, Noriyuki Ohi, Tatsuo Nishiki, Hikaru Inoue, Kazuo Minami, Seiya Nishizawa,
Miwako Tsuji, Naonori Ueda, 152K-computer-node parallel scalable implicit solver for dynamic nonlinear
earthquake simulation, HPCAsia2022: International Conference on High Performance Computing in Asia
Pacific Region, 18-29, 2022.

2022/10/27

The state-of-the-art solvers@Supercomputing conference

Decomposed to two
independent small problems
Interactions are cut off…

+

Unprecedented & 10 exa-scale problem

一般化可能な新たなアルゴリズム:
Integration of HPC & Data-learning is developed!
Equation-based approach + Data-driven approachwith High performance computing on Supercomputer+
京から富岳への42倍のハードウェア性能向上

10 exa-scale の解析性能, Ground-breaking solver!

2022/10/27

Equation-based と data-driven approach• Equation-based approach• 低次モードは低コストで構成可能だが，高次モードの構成は高コスト• Data-driven approach• コストと精度は用いるデータサイズによる．
• 精度保証は一般に難しい• 全領域の高次モードの構成は難しいが，小領域ではある程度までは可• Equation-based と data-driven approachを組み合わせることで効率的に広帯域の解を構成• Equation-based は低次モードの構成に，data-driven approachは高次モードの構成に用いて，これを初期解推定とする

Equation-based + data-driven initial solution estimator• In dynamic problem• Estimate initial solution → solve Ax=f with iterative solver• Resolves high-order modes prior to iterative solver using resultsof previous time-steps• Time evolution matrix 𝐴 is estimated as 𝑋
= 𝐴𝑋
, where 𝑋
=
𝑥
,𝑥
,…,𝑥
using previous 𝑠 steps of 𝑥 = 𝛿𝑢 − 𝛿𝑢• Transient low-order modes that deteriorates accuracy of data-driven
method is removed using equation-based predictor (𝛿𝑢by AdamsBashforth method), and data-driven approach is used to obtain nontransient high-order modes:
𝛿𝑢 = 𝛿𝑢+ 𝐴 𝛿𝑢− 𝛿𝑢• Highly scalable on large systems as computation is local without
communication• Time evolution matrix 𝐴 can be constructed/evaluated per subdomaininternal to MPI process domains

2022/10/27

Multi-grid preconditioner• Most high-order modes resolved with data-driven initial solution predictor;
remaining high-order mode errors concentrated in parts with especially poor
convergence• Design W-cycle multi-grid preconditioner that specifically refines high-ordermodes at parts with poor convergency
Solve preconditioning matrix equation 𝐴𝑧 = 𝑟 roughly(single precision)
CGprehmr
CGpre
CGpre

CG loop

Second ordered
tetrahedron

Solvingpreconditioning
matrixComputations ofouter loop

Equation to be solved(double precision)

Part ofsecondorderedtetrahedron
mesh

Second ordered
tetrahedron 𝐴𝑧 = 𝑟

Outer loop
CGprelmr• This new algorithm is applied toearthquake simulation• The 1,070-fold capability
improvement (25.45(=9318/366) X
42) far beyond the state of the art
is achieved!

• Low-order unstructured nonlinear
dynamic finite element problem(1.201 trillion DOF) solved by implicit
solver
• 19.216 trillion component dataset
generated inside simulation arelearned every time-step for improving
simulation performance

Linear tetrahedron
𝐴 𝑧 =𝑟CGprelmr

State of the art
Failed in this problem

developed※全て富岳上での計測

2022/10/27• Three-dimensional low-order unstructured nonlinear
finite element analysis of an inland earthquake on a
fault-structure model comprising a highly integrated
city near a terminal station with soil and crust of Kanto
Plain in Japan.

10 exa-scale capability• The application is an ultra-large-scale and ultra-highfidelity problem: the domain size of the entire model is
256000 × 205000 ×100000 m, discretization size
varies from approximately 0.125 to 64 m, the degree of
freedom is 324,006,449,076, and the number of
second-order tetrahedral elements is 80,540,873,752.
Enlarged view shows responses of the city with
underground structures.

• Characteristics of target matrix is not so good…

• analyze up to 9550 time-steps (dt=0.005 sec) on Fugaku• With full system (152,352 computer nodes (609,408 MPI processes × 12 OpenMP threads =7,312,896 parallel computation)) in 0.77 sec per time step

まとめ
• 計算科学• 物理シミュレーションの性能向上を目指した従来型のシミュレーショ
ンに加え，データとコンピューティングの連成を陽に考慮することで，解析能力を向上しうる可能性• 演算の質・計算アーキテクチャの性質を考慮することで，従来とは異なるタイプのアルゴリズム開発の可能性• 地震シミュレーション• 従来難しいとされてきたシミュレーションがいくつか実現され，それに応じたサイエンス的検討の可能性• サイエンス側からの新たな計算科学への要請