Title pattern 1

(1)

都市の地震災害シミュレーションのための

非線形波動場解析手法の開発

理化学研究所計算科学研究機構

(

日本学術振興会

PD)

藤田航平

2014

年

12 月

19 日

(

金

)

平成

26 年度「京」における高速化ワークショップ

(2)

発表の概要

•

1. 地震シミュレーションにおける，非線形・非構造有限要素法の高速化の必要性

2. 高速化方法

(

アルゴリズムを中心に

)

3. 高速化によってできるようになること

• 共同研究者

• 東京大学地震研究所市村強，田中聖三，堀宗朗，

Lalith Maddegedara

• 高度情報科学技術研究機構志沢由久，小林寛

• 詳細は

”Physics-based urban earthquake simulation enhanced by 10.7 BlnDOF

x 30 K time-step unstructured FE non-linear seismic wave simulation,

(3)

はじめに

• 地震災害

• 物理は単純

• ただし，領域サイズ・分解能が高く，非線形性が強い

• 物理を

3 次元で追うことは大規模な解析に帰着する

• 従来は統計的・経験的手法で対処

• 物理を忠実に解くには，

HPC

による大規模解析が必要

過去の地震被害データに基づく統計・経験式震源からの距離揺れの強さ揺れの強さ構造物の被害確率地震災害の物理 – 波動方程式

 

2₂ 2 2 x u u c t u     

(4)

地震シミュレーションの現状

地殻解析 地盤解析 構造物解析 Resolution O(0.1km), Domain size O(100km) Resolution O(1m), Domain size O(1km) Resolution O(0.1m), Domain size O(10m)

•

3 段階の組み合わせ

• HPC

アプリケーションでは「地殻」と「構造物」が主な対象

• 「地盤」の解析が重要かつチャレンジングな問題に

– 周波数成分によって地震動に10倍の増幅がかかる – 構造物と人命の安全に直接かかわる – 地震シミュレーション全体の信頼性向上に不可欠

(5)

「地殻解析」と「地盤解析」の違い

• 地盤解析では地殻解析と異なる問題を解く必要がある

• 性能は出しにくいが，重要な問題

「地盤」中の波動伝播

• 複雑形状 • 多数の層からなる • 非常に柔らかい土からなる • 非線形の波動方程式を解く • 安定条件が厳しい – 陰的な時間積分が適している – 方程式の求解＋集団通信が必要 O(10m)

 

2₂ 2 2 x u u c t u     

「地殻」中の波動伝播

• 形状はそこまで複雑ではない • 少数の層からなる • 固い岩からなる • 線形の波動方程式で近似 • 安定条件はそこまで厳しくない – 陽的な時間積分が広く使われている – 方程式の求解は不要，集団通信も不要 2 2 2 2 x u c t u     

(6)

地盤解析の現状

[1] Ichimura, T. et al., Journal of Pressure Vessel Technology, 2014.

• 陰的時間積分・低次要素を使った有限要素法による地盤解析

[1]

– 空間分解能

4m,

時間分解能

2.5Hz

• 計算負荷が大きい

– 構造物の固有周期である

10

1

_Hz

_{までは計算できない}

→

_{構造物解析の入力として使うこと}

ができない

– 多数ケースの解析を実行できない

→

解析結果の信頼性確認が難しい

(7)

ターゲット問題

• 大自由度・陰解法・非構造格子有限要素法の高速化が地盤解

析の実用化に必須

➥

時間分解能

15Hz

で実問題を解くには，

マトリクス

(

毎タイムステップ変化する

)

未知ベクトル(100億自由度)

外力ベクトル

30,000

ステップある時間ステップ毎に相対誤差

10

-8

_で求解

(

各タイムステップで物性を線形化して解析

)

(8)

アルゴリズムによる高速化

• マルチグリッド法

– 計算量を削減

• 精度混合演算

– 計算コストを削減

– 要求

B/F

を減らす

• Etc…

• これらを組み合わせてより高速解法を開発

http://www.mgnet.org/mgnet/tutorials/xwb/mg.html

マルチグリッド法

(9)

GAMERA

• 非線形地盤解析用のソルバー

• HPC

手法を組み合わせることで開発

–

multi-

G

rid method

– A

daptive conjugate gradient method

– M

ulti-precision arithmetic

– E

lement-by-element method

–

p

R

edictor with

A

dams-bashworth

method

(10)

GAMERA

のアルゴリズム

• 前処理行列を荒く・高速に解くことで，

CG

ソルバーのループ数を削減

– 前処理にマルチグリッドを使用することで，前処理の計算コストを削減 – 前処理に単精度を使うことで計算と通信を減らす

Equation to be solved

(second ordered tet. mesh, double precision) CG loop

Computation of outer CG loop

Outer loop

Solve system roughly using CG solver (linear tetrahedral mesh, single precision)

Solve system roughly using CG solver

(second ordered tet. mesh, single precision) Use as initial solution

Use for preconditioner of outer loop Solving preconditioning matrix Inner coarse loop

Inner fine loop

Solving

preconditioning matrix

(11)

GAMERA

のアルゴリズム

•

1 ループあたりの通信量

(byte

比

)

Outer

: Inner_fine :

Inner_coarse

=

1 : 1/2 :

1/8

•

1 ループあたりの計算量

(FLOP

比

)

Outer

: Inner_fine :

Inner_coarse

=

1 : 1 :

0.18

• Inner_fine

と

Inner_coarse

は単精度計算

→

要求

B/F

が低く，キャッシュを実質的に

2 倍使える

Equation to be solved

(second ordered tet. mesh, double precision) CG loop

Computation of outer CG loop

Outer loop

Solve system roughly using CG solver (linear tetrahedral mesh, single precision)

Solve system roughly using CG solver

(second ordered tet. mesh, single precision) Use as initial solution

Use for preconditioner of outer loop Solving preconditioning matrix Inner coarse loop

Inner fine loop

Solving

preconditioning matrix

(12)

計算環境

• 京コンピュータ＠理化学研究所・計算科学研究機構

– 82,944 compute nodes, each with single SPARC64 VIIIfx CPU, connected with Tofu interconnect (6-dimensional mesh/torus network)

– Peak performance of each CPU: (FMA SIMDs with vector length 2) x 2 sets x 8 cores x 2.0GHz = 128GFLOPS

(13)

性能測定

• 問題設定

– x

,

y方向に並べることで，周期性のある問題を設定

–

METIS

で領域分割：実問題と同様のロードバランス特性

Layer 1 2

Pressure wave velocity (m/s) 700 2,100 Shear wave velocity (m/s) 100 700 Density (kg/m3₎ _1,500 _2,100 Damping 0.25 (h_max) 0.05 Strain Criteria 0.007 - Material properties 60m 20m 64m 8m Layer 1 Layer 2 x y z Element size: 1.0m

(14)

• 京全系を使ったテスト

(82,944

ノード

, 27 billion DOF)

– PCG法と比べて，「反復一回当たりの計算時間短縮」＋「トータルの反復回数が少なくなる」ことで， 3.53倍の高速化

アルゴリズムの効果

0 1 2 3 4 5 6 7 Inner coarse Inner fine Outer Elapsed time for solving 1 time step (s)

(125) (68) (9) (125) (68) (9) (8) (399) (355) () ：収束までに要した反復回数 *3x3のブロックヤコビスケーリングを前処理に使用したCGソルバー GAMERA Without mixed precision Without mixed precision and multigrid PCG method* Elapsed time (s) 1.63 1.87 6.95 5.75 FLOPS/Peak (%) 7.57 6.60 7.88 8.11

(15)

GAMERA

の高速化

• ここまでのまとめ

– 通信・計算を減らすアルゴリズムを使うことで，京全

系を使った問題において従来手法に比べて

_3.5

倍高

速化

• ここからの概要

– 行列ベクトル積の計算方法の変更と，スレッドチュー

ニングでさらなる高速化を実施

(16)

行列ベクトル積のアルゴリズム変更による高速化

• 通常の計算方法 (Compressed Row Storage法，等)

• マトリクスは大きすぎてキャッシュにのらず，メモリアクセスが多い

• 非線形性によりタイムステップ毎にCRS形式でマトリクスを準備しなおす必要あり

• Element-by-Element (EBE)法による計算

• 行列ベクトル積を計算するたびに要素マトリクスを計算

• マトリクスを記憶する必要がなく，右辺・左辺ベクトルは小さいので，キャッシュにのせることができる

• 以前はouterとinner_fineはEBE法，inner_coarseはCRS法だったところを，すべてEBE法に変更

• 計算量は増えるが，「ピーク性能比の向上」＋「CRSの準備部分が不要」となり，全体として高速化行列ベクトル積 EBE法 (各要素に対して) = … += … CRS法 (各自由度に対して) = ₊₌ 要素マトリクス =

(17)

EBE

法のスレッド足し込みの改善

変更前

Allocate temporary vectors per

OpenMP thread, and initialize to 0 Update components by EBE (red) Add all components

変更後

Initialize components necessary Update components by EBE (red) Add components necessary + + = + + = + + = + + = + + = + + = + = = = + = = + =

• Element-by-Element

法による左辺ベクトルの足しこみ

• 図：ベクトルの長さ

6,

スレッド数

3 の場合

(18)

0 50 100 150 高速化前高速化後 others inner coarse inner fine outer

変更の効果

計算時間(全体，s) 158.9 99.6 計算時間(ソルバー部，s) 129.6 91.5 ソルバー部の実行効率(FLOPS/Peak, %) 9.00 14.49 ソルバー部のメモリスループット(%) 13.23 6.89

Elapsed time for solving 100 steps (s)

• Inner_coarse

部：「

CRS

準備

+ CRS

計算」より，

EBE

計算の方が速くなった

• Inner_fine, outer

部：スレッド足しこみの改善などにより高速化

• 全体として

158.9s → 99.6s (1.59

倍の高速化

)

EBE計算 CRS準備+ CRS計算京4,608ノードで計測全体ソルバー

(19)

ストロング・スケーリング

• 問題規模

(120

億自由度

)

を一定に保ったまま計算ノード数を拡大

CPUコア数 (計算ノード数) 36,864 (4,608) 73,728 (9,216) 147,456 (18,432) 294,912 (36,864) 計算ノードあたりの自由度数 2612k 1306k 653k 327k ソルバー部の経過時間 (s) 686.5 306.3 180.0 111.0 ソルバー部の実行性能 (ピーク比%) 15.4 17.1 14.7 11.8

Elapsed time for solving 100 steps (s) # of CPU cores 64 256 1024 Total Solver

(20)

ウィーク・スケーリング

• 計算ノードあたりの問題規模を一定に保ったまま計算ノード数を拡大

0 40 80 120 Others Inner coarse Inner fine Outer CPUコア数 (ノード数) 36,864 (4,608) 73,728 (9,216) 147,456 (18,432) 294,912 (36,864) 663,552 (82,944) 自由度数 15億 30億 60億 120億 270億 CPUコアあたりの自由度数 41k 41k 41k 41k 41k ソルバー部の実行性能 (ピーク比%) 14.49 14.06 13.28 11.82 12.34

Elapsed time for solving 100 steps (s) T ot al elaps ed tim e Solv er

(21)

適用例

• 新宿を対象に，

1995

年

JMA

神戸で観測された波形を入力

Acceleration NS

Acceleration EW

Acceleration UD

(22)

高速・大規模解析でできること

• シミュレーション結果の信頼性確認

(Verification and Validation)

–

Verification:

物理モデルが正確に解かれているか，

Validation:

シミュレーション結果が実際の現象に一致しているか，

に分けてシミュレーション結果の確からしさを検証する

物理現象物理モデル(支配方程式) シミュレーション結果モデル化数値解析 Verification Validation

An Overview of the PTC 60 / V&V 10, Guide for Verification and Validation in Computational Solid Mechanics, ASME, http://cstools.asme.org/csconnect/pdf/CommitteeFiles/24816.pdf

(23)

• 従来は，解析領域の一部を取り出した準備解析により，メッシュサイズを

決めていた

– 局所的に揺れが増幅する箇所では収束が甘くなる可能性

• 同じ数理問題に対して，分解能を変えた解析モデルを作成し，解析領域

全体の解の収束性を確かめる

高速・大規模解析でできること

対象周波数

要素サイズ

自由度数

15Hz

0.66m

107 億

10Hz

1.0m

33 億

5Hz

2.0m

4.7 億

(24)

解の収束性確認

• 15Hz

モデルを使うことで，地表面での地震動分布は

1%

以内の誤差で収束

15Hz model 15Hz – 5Hz model 15Hz – 10Hz model SI (m/s) Difference of SI (m/s)

(25)

解の収束性確認

• 時刻歴波形も

1%

以内の誤差で収束

• 京＋

GAMERA

によって初めて都市スケールで実現

Convergence of structure response Convergence of ground acceleration

A cce lerat ion a t su rf a ce (E W , in g a l) Di sp lace m e n t a t roo f (E W , in m ) t (s) 15Hz model 15Hz – 10Hz 15Hz – 5Hz 15Hz model 15Hz – 10Hz 15Hz – 5Hz

(26)

Title pattern 1

都市の地震災害シミュレーションのための

非線形波動場解析手法の開発

理化学研究所 計算科学研究機構

(

日本学術振興会

PD)

藤田 航平

2014

年

12

月

19

日

(

金

)

平成

26

年度「京」における高速化ワークショップ

発表の概要

•

目次

1.

地震シミュレーションにおける，非線形・非構造有限要素法の高速化の必要性

2.

高速化方法

(

アルゴリズムを中心に

)

3.

高速化によってできるようになること

•

共同研究者

•

東京大学地震研究所 市村強，田中聖三，堀宗朗，

Lalith Maddegedara

•

高度情報科学技術研究機構 志沢由久，小林寛

•

詳細は

”Physics-based urban earthquake simulation enhanced by 10.7 BlnDOF

x 30 K time-step unstructured FE non-linear seismic wave simulation,

はじめに

• 地震災害

•

物理は単純

•

ただし，領域サイズ・分解能が高く，非線形性が強い

• 物理を

3

次元で追うことは大規模な解析に帰着する

•

従来は統計的・経験的手法で対処

•

物理を忠実に解くには，

HPC

による大規模解析が必要

 

地震シミュレーションの現状

•

3

段階の組み合わせ

•

HPC

アプリケーションでは「地殻」と「構造物」が主な対象

•

「地盤」の解析が重要かつチャレンジングな問題に

「地殻解析」と「地盤解析」の違い

•

地盤解析では地殻解析と異なる問題を解く必要がある

•

性能は出しにくいが，重要な問題

「地盤」中の波動伝播

 

「地殻」中の波動伝播

地盤解析の現状

•

陰的時間積分・低次要素を使った有限要素法による地盤解析

– 空間分解能

理化学研究所計算科学研究機構

藤田航平

東京大学地震研究所市村強，田中聖三，堀宗朗，

高度情報科学技術研究機構志沢由久，小林寛

_Hz

_{までは計算できない}

_{構造物解析の入力として使うこと}

_で求解