災害軽減へ向けて
理学・工学をつなぐ計算科学
市村強
東京大学地震研究所
計算地球科学研究センター

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

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

数理問題の観点からみると，見通しがつきやすい（レングスケー
ルや分解能が違う位で，理と工をそこまで区別する必要はない）
・これらの順・逆解析をする問題に帰着
・問題によって「構成則」や「感覚」は違うが．

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

日本列島規模784km

y

町内会？規模1184kmzx

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

東日本規模

400km

地下高速道路トンネル

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

理学・工学をつないで災害軽減へ
観測

解析

現象の理解

投資
災害軽減

知見の利活用

目的（応力フリーの）境界条件や幾何形状が解に強く影響を及ぼす数理問題なので，有限要素法を使いたい．これが出来れば，理も工も共通の解析基盤にのる？
しかし，対象領域が大きく，必要分解能が高いため，解の収束性（信頼性）を担保しようとすると大規模問題（自由度108-12）になってしまい，解けない．結果的にsimplified モデル（半無限近似，水平成層，近似数値解法＋物理的な近似の解釈＝他分野に分かりずらい）を用いざるを得なかった．
様々な問題の解析基盤となる大規模な有限要素解析を実施可能なように（計算科学的検討を精力的に行うことで）解析高度化をめざす
ちなみに，固体を扱う製造業も「有限要素法で計算する」ことになり，かなりの投資がなされ，有限要素法のソフトが整備されてきた．ただし，工業製品くらいの大きさなので，その対象となる自由度は100万自由度程度．

100万並列を可能とする超並列計算用ソルバー
Solving preconditioning
matrix
Inner coarse loop

CG loop
Solvingpreconditioning
matrixSecond ordered
tetrahedron

Computations ofouter loop

Equation to be solved(double precision)

Outer loop

Solve system roughly using CG solver(single precision)

Inner fine loop

Use as initial solution

Linear tetrahedron

Solve system roughly using CG solver(single precision)
Use for preconditioner of outer loop

Second ordered
tetrahedron• Solve preconditioning matrix roughly to reduce number of CG loops• Use multi-grid method to reduce cost of preconditioner• Use single precision in preconditioner to reduce computation & communication

Resource to solve target problem• K computer @ RIKEN Advanced Institute for ComputationalScience• 10.6 PFLOP machine consisting of 8 cores x 82944 compute nodes• Maximum available time: 8.0 hours• Godzilla can solve target problem using full machine within
available time!

• Attain superb scalability for utilizing all 663,552 cores• Attain high efficiency & performance for finishing simulation in 8.0 hours

Performance of this Solver (1.08 trillion DOF @full K computer)• 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) on full K computer• 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

2
Normalized
9216
elapsed time 1
1/2
1/4
1/8
1/16
1/32

294912 663552# of coresSize-up scalability
20

Weak Scaling on Summit @ SC18 solver
500
400

349.8

374.6

380.2

415.1

454.2

311.7

327.3

75.8

77.680.4

82.9

84.3

83.7

90100.4# of GPUs

288
288

576
576

1152
1152

2304
2304

4608
4608

6144
6144

12288
12288

24576
24576

# of node

48

96192

384

768

1024

2048

4196

DOF

3.5×109

7×109

14×10928×10956×10975×109151×109

302×109

SC18’s efficiency toFP64 peak

19.5%

SC18MOTHRASC14GAMERA

Elapsed time (s)

302.5

300
200
1000

14.7%

Difficulties in FEM for huge problems• Analysis Cost

• Solver on supercomputer• Scalability communication between compute nodes for synchronize values(vector-vector, matrix-vector product)• Cost for Constructing FE model• Robust & good aspect ratio & fast

2 types of Codes for Construction of FE model• For layer media• Up to Trillion-order(1012-) DOF• Unstructured
tetrahedral elements• with structured grid
mapping• For more complicated• Up to billion-order(109-) DOF• Unstructured
tetrahedral elements• With multi-material
marching cubic

地盤・プラント丸ごとシミュレーション
要素サイズ16㎜の四面体二次要素
自由度約490億24

余効変動データを用いた粘弾性構造・余効すべりの同時推定手法：
数十億自由度の高詳細有限要素モデルを用いたパラメータ最適化
η1
η2

η3

東北地方の高詳細有限要素モデルとその断面
粘性層真値(Pa s)

推定値(Pa s)

η15.000×1018

5.012×1018

η25.000×1018

4.052×1018

η31.000×1018

1.020×1018

粘弾性真値と推定結果• 数十億自由度の高詳細有限要素モデ
ルを用い，アセノスフェア粘性率と余効すべりの同時推定手法開発• 高速ソルバーとアジョイント法を組み合わせた計算コスト的に効率の良い非線形最適化• 東北地震を模擬した数値実験にて良好な結果データとも整合する，より信頼性の高いシミュレーションが可能になる？

真のすべり量（左：地震時，右: 余効すべり）

推定すべり量プレート形状の曖昧さを考慮したインバージョン：8000万自由度の地殻変動解析を36万回実行• 高詳細な8000万自由度の有限要素地殻モデルを1000種類構築し，小断層180＊2成分に対して地殻変動解析を行いグリーン関数を評価• 8000万自由度の地殻変動解
析一回計算するのに，地震研の小規模GPUクラスタで20秒しかかからなくなった
FE model and locations of grid points used as inputs for the basis functions.
• 1000種類の8000万自由度の有限要素モデルを自動生成

曖昧さを陽に考慮した評価が可能になる？

地殻応力詳細解析のための超高分解能地殻変動解析：プレート面上で101-2mの分解能（２兆自由度）
Target domain

Macro-analysis region

Input slip distribution(one case of data assimilation simulation)

Micro-analysis region

Compute displacement and stress distribution• 列島全域の高詳細地殻モデ
ルを階層型解析により分割• 列島規模の100億自由度有限要素地殻モデルと局所的な2兆自由度有限要素地殻モデ
ルを生成し，地殻変動解析を実施．

Continentalplate

プレート面上での応力評価が可能になる？

Pacificplate
Philippine
sea plate

Macro-analysis region (0.01 T-DOF
problem with 500 m resolution in [1])Micro-analysis region (2.05 T-DOF
problem with 45 m resolution)

First ever city simulation using full power of K computer• Target Tokyo metropolis• One of world’s largest
metropolitan areas• Enabling of this area leads toapplication to nearly all cities

5 miles (8.04km)

Fault-to-city simulation

29

Fault-to-city simulation• Earthquake linear wave propagation simulation @ 10Hz target frequency• Use 18 billion unstructured low ordered finite elements & 56 billion DOF model with
implicit time integration

UNEFault

City simulation
domainVelocity Norm (cm/s)
10.0
5.0

T=15 s

T=16 s

T=17 s0.030

City simulation

31

PDE求解と情報科学の融合：人工知能により系を学習し，方程式求解を高速化• 微分方程式の持っている情報を用いる• 方程式の解きにくい部分の情報を支配方程式がもっている• 解きにくい部分をAIで抽出して，優先して解いてしまう

Phenomena

Dataanalytics

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

Ax=f• 数理問題という観点からは理工連携の解析は見通しがつきやすい．それを下支えする大規模有限要素法も計算科学的検討を精力的に行うことで大幅に進歩した．

まとめ• 理工の共通の解析基盤が出来れば，左図のサイクルがさらによくまわると期待．

観測• 解析基盤だけでなく，データ解析にも計算科学は有効．うまくいけばさらによくまわるようになる
現象の理解と期待．

解析

投資

知見の利活用

災害軽減
デジタルツイン地球多数の観測機器から得られたデータを蓄積し，大規模解析基盤を用いて多種多様なシミュレーションを駆使して過去・現在・未来を評
価可能なデジタルツインを観測データとの融合により構築．理学から工学をシームレスに接続し，地震・火山に関する諸現象の科学的研究と直接又は間接に地震・火山に起因する災害の予防並に軽減方策を探究する• 解析の高度化により，観測の意義がより明確化されるとともに，理工連携がさらにすすみ災害軽減への寄与へ• なお．．．
• supercomputerを使えば，巨大な計算が簡単にできるわけではなく，超並列技術が必須．
• 1億自由度くらいであればサーバーくらいで高速計算できるバージョンもあります．
• 今回は固体を中心に構成．固体以外も同様の傾向