2014年10月29日
サイエンティフィック・システム研究会
科学技術計算分科会 2014年度会合
次世代HPCを支える技術
ホーテルオークラ神戸

京コンピュータを使った地盤震動
シミュレーションの現状
理化学研究所 計算科学研究機構(日本学術振興会PD)
藤田航平

1

地震・津波災害： 現代日本における主要な自然災害~1960：地震・津波＋風水害

1960~：地震・津波

http://www.cas.go.jp/jp/seisaku/resilience/dai1/siryou6-1.pdf••

自然災害による人的被害の多くは防げるようになったが，地震・津波による被害は50年間解決できていない
初期のスパコンの誕生も1960年ごろ．この50年で飛躍的に発展したHPCで地震・津波問題を解決できないか？
2

はじめに (1)• 計算科学の観点から見た地震・津波災害– 物理現象自体は単純(e.g., 波動方程式)– ただ，領域が広大な割に必要な時空間分解能が高い• 従来の地震被害予測– 物理を3次元で追うことは膨大な計算となるため，統計手法や経験式を組み合わせて予測– 観測データが少ないため大数の法則が効きにくく，信頼性は必ずしも高くない
揺れの強さ
距離減衰式

震源からの距離
構造物被害の確率

地震の物理は，連続体中の波動伝播

フラジリティ・カーブ

揺れの強さ
3

はじめに (2)
• 物理を3次元で解くことで，地震災害の現象を理解や予測の信頼性を高めることができると期待– HPCのアプリケーションとして世界中で開発が進められている• HPCI戦略プログラム分野3での「地震・津波の予測精度の高度化に関する研究」– 広域の問題を大規模計算で正確に解く– 一連の災害現象をつなげて解く (統合地震シミュレーション)➥ 京の能力・HPC技術を実社会の安全に反映するためのミッション
City

Structural response

Surface soil

Bedrock
Crust
Fault

Soil amplification
Wave propagation4地震シミュレーションの現状• 地震災害の3段階を組み合わせて解析• 従来から「地殻中の波動伝播」と「構造物解析」はHPCアプリケーションとして盛んに研究されてきた• 都市スケールの地盤震動問題が重要かつチャレンジングな問題となっている– 構造物の地震被害に大きく影響• 特定の周波数成分によっては10倍程度の増幅がかかることも– 広域・複雑形状・非線形物性のため，特に大規模な問題になる– 地震シミュレーション全体の信頼性を上げるために不可欠な要素
分解能O(0.1km), 領域O(100km)

分解能O(0.1m)
領域O(10m)
分解能O(1m), 領域O(1km)

地殻中の波動伝播[Furumura et al.]

都市スケールの地盤震動

構造物解析
5

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

地盤震動問題– 地表～深さ100m程度の地盤は，堆積や浸食などで複雑な物性分布・形状となっている– これを正確に解くには，大規模unstructured FEM解析が適している•

都市スケールの3D地盤震動シミュレーションは可能になりつつある– 空間分解能4m, 時間分解能2.5Hz•

その一方で，計算コストが大きいため，– 構造物の応答に支配的な分解能(15Hz)は達成できていない– 複数回の計算ができないため，結果の信頼性を確かめることも難しい•

高速な大規模unstructured FEM解析が地盤震動シミュレーションを実用につなげるために不可欠[1] Ichimura, T. et al., Journal of Pressure Vessel Technology, 2014.

6

大規模unstructured FEM解析の必要性• 地盤震動問題だけでなく，地震・津波問題全般において，「形状を正確に解く」ことが重要になっている– 津波入力の推定精度の向上：形状を正確に反映した計算で，観測データを使った断層状態の推定精度向上が期待できる• 「大規模unstructured FEM解析」は地震・津波問題の基盤的な手法に– 陰解法・非構造格子という性能の出しにくい問題だが，意義の大きい問題

DONET, DONET2: 地震・津波観測監視システム
https://www.jamstec.go.jp/donet/e/donet/donet2.html

7

大規模unstructured FEM解析• 陰的な時間積分を使うため，疎行列線形方程式ソルバーが重要• 機械・原子力などの工学分野で長く研究されてきた– シェル・パイプなどの薄い・細い部材などにより，性質の悪い大規模問題となる– FrontISTR, ADVENTUREなど高性能並列解析ソフトウェアが開発• 地震・津波災害系の問題：– 地盤・地殻といったバルクのある問題のため，性質はそこまで悪くないが，より大規模な問題– 本研究はここを狙う

車体：大部分はシェル部材の組み合わせ

地殻・地盤：中の詰まった固体
8

問題設定• 離散化された波動方程式を解く– 非線形物性のため，各ステップでマトリクスが変化– 毎ステップ，大規模な対称疎行列の線形方程式を解く問題に帰着• マトリクスが変化するため，高性能前処理はコストが見合わない• その一方で，静的問題よりも対角項が大きくなるので，マトリクスの性質はそこまで悪くない
Unknown vector

Matrix (componentschanges every time step)

Force vector
9使えそうな方法• マルチグリッド法– 計算量を減らす• 精度混合演算– 計算コスト・要求B/Fを減らす• Element-by-element法– 要求B/Fを減らす• など…

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

マルチグリッド法• それぞれ有用だが，組み合わせることでこの問題をより高性能で解けないか
10

使えそうな方法• マルチグリッド法– 計算量を減らす• 精度混合演算– 計算コスト・要求B/Fを減らす• Element-by-element法– 要求B/Fを減らす• など…

b  AuMatrix-vector productb   A e x e , c ,  uee(for each element)

Matrix-vector productcomputed by elementby-element method

Element-by-element法:
要素毎に行列ベクトル積を計算・足込み• それぞれ有用だが，組み合わせることでこの問題をより高性能で解けないか
11

GAMERA• これらのHPC手法を組み合わせて，地盤震動問題を高速に解く– multi-Grid method– Adaptive conjugate gradientmethod– Multi-precision arithmetic– Element-by-element method– pRedictor with Adams-bashworth
methodTsuyoshi 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, SC’14, Gordon Bell Finalist, in press.

12

GAMERA, solver for wave amplification problem• Solve preconditioning matrix roughly to reduce number of CG loops– Use multigrid method to reduce compute cost of preconditioner– Use single precision in preconditioner to reduce computation and communication• Parallelized using MPI, OpenMP, auto parallelization
Solving preconditioning matrix
CG loop
Solvingpreconditioningmatrix
Othercomputationsof CG loop
Outer loop
Equation to be solved(second ordered tet. mesh,double precision)

Inner coarse loop
Solve system roughly using CG
solver (linear tetrahedral mesh,single precision)
Inner fine loop

Use as initial solution

Solve system roughly using CG
solver (second ordered tet. mesh,single precision)
Use for preconditioner of outer loop

13

Performance measurement• Problem setting– 94,407,951 DOF with element size of 1 m– Duplicated in x, y directions to make problems with periodical geometry

14

Size-up scalability of GAMERA• 計算ノードあたりの問題規模を一定に保ったまま，並列数を増やす– 非線形波動場解析を100ステップ解くのにかかった時間を計測– 陰解法・低次要素・非構造格子という性能の出しにくい問題だが，京全系でsolveできている

120

OthersSolverTotal elapsedtime

Elapsed time forsolving 100 steps (s)Inner coarse
Inner fine80
40Outer0

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コアあたりのDOF数

41k

41k

41k

41k

41k

FLOPS/Peak of solver (%)

14.49

14.06

13.28

11.82

12.34
15

Speed-up scalability of GAMERA
120億DOFの問題を並列数を変えて計測
512
Elapsed time (s)

256

Total

128

Outer

64

Inner fine32

Inner coarse

16
8

73,728(9,216)

147,456(18,432)

294,912(36,864)

CPUコアあたりのDOF数

163k

81k

41k

Elapsed time of solver (s)

306.3

180.0

111.0

FLOPS/Peak (%)

17.1

14.7

11.8

CPUコア数(ノード数)

Elapsed time (s)

128
64

32
16842TotalComputation
Communication
Synchronization
16Effectiveness of algorithm• GAMERAのアルゴリズムの効果を計測– 京全系を使って計測 (82,944ノード，270億DOF)
Elapsed time (s)() indicates the number of loops
elapsed for convergence400300Inner coarse
Inner fineOuter

(28864)200
1000(25391)(10065)(6691)( 614)GAMERA(10166)(6658)( 619)( 609)

Without
mixed
precision

Withoutmixed
precision andmultigrid

PCG
method*Elapsed time (s)

113.2

131.8

344.4

314.4

FLOPS/Peak (%)

11.7

10.2

10.8

10.0

Mem. Throughput/Peak (%)

5.59

10.3

12.0

12.4*Use 3x3 Block Jacobi method for preconditioning

17

Portability of GAMERA• Intel系のCPUを使ったStampede を使ってGAMERAのポータビリティを検証– 単精度計算の速い計算機を使うことで，精度混合演算がより効果的になる

18

適用例：問題設定# of stories• 3層の地盤構造・13,275棟の構造物からなる2.0 x 2.0 kmの領域(新宿付近)• 地盤工学会提供のデータをもとに地盤モデルを作成

LayerStructures

Surface topology

123

Vp (m/s)

1,210

1,380

1,770

Vs (m/s)

150

255

490

1,500

1,800

1,900

0.25 (hmax )

0.05

0.005

0.005

--

Density (kg/m3)
Damping

Elevation (m)Strain Criteria

Interface between layer 1 & 2Bedrock19

適用例：解析モデル• 要素サイズ0.66 mでモデル化することで，10,755,536,091 DOF• 神戸波を入力 (dt = 0.001 s, 30,000ステップ)

Structure A

Structure BNUE20

計算コスト• 京コンピュータ36,864ノードを使って，7時間44分 (27,837s)– 91%の時間がソルバーに使われており，実問題でも高性能な計算が実行できている
Elapsed time for 30,000 steps (s)Others2,539.0
Others:- Compute nonlinear constitutivemodel- Compute block diagonal matrix- Read input boundary condition- Write output data

Solver25,298.0

Total: 41,521 s
Sustained FLOPS/peak: 12.6% (solver), 11.8% (Total)
21

Magnitude of grounddisplacement (m)00.3

Magnitude of relative
displacement of structures (m)
00.5

解析結果

00.7
Maximum relative
300
100
displacement of structures (m) Horizontal magnitudeof SI at surface (cm/s)

t = 10.0 (s)t = 9.5(s)t = 9.0 (s)Time history response of ground and structures

GroundmotionStructure
response• 地盤と構造物の特性を踏まえた被害分布を求めることができる
22

シミュレーション結果の信頼性確認(1)• V&V (Verification and Validation)[3]– Methodology to ensure the reliability of simulation results– Consists of two phases• Verification: check if the physics model is solved correctly• Validation: check if simulation results follow the actual phenomenon
Physical phenomenon
Physical modeling
Validation

Physics model (governing equation)
Numerical analysis

Verification

Simulation results[3] 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

シミュレーション結果の信頼性確認(2)• 従来は，解析領域の一部を取り出した準備解析により，メッシュサイズを決めていた– 非線形な挙動で局所的に揺れが増幅する箇所では収束が甘くなる可能性がある• 同じ数理問題に対して，分解能を変えた解析モデルを作成し，解析領域全体の解の収束性を確かめる– 5Hz (ds=2m), 10Hz (ds=1m), 15Hz (ds=0.66m)の3つのモデル– シミュレーションの信頼性向上に重要なステップ

24

解の収束性確認• 地表面でのSI値の分布の収束を確認– SI値：構造物の被害を見積もるのに使われる地震動強さの指標– 15Hzモデルで，誤差1％以内に収束m/s
15Hz model

15Hz – 10Hzm/s
15Hz – 10Hz25

解の収束性確認• 加速度波形も15Hzモデルで誤差1%以内に収束
Point A

800gal15Hz model
30gal
30gal

Point B

15Hz - 10Hz
15Hz - 5Hz

1250gal15Hz modelm/s
15Hz model のSI値分布

25gal

15Hz - 10Hz25gal

15Hz - 5Hz

0t (s)
30
Acceleration waveform at surface (East-West direction)26

構造物応答の収束性確認

Structure A displacement EW (m)• 2階建てRC構造物2棟の変位を比較– 5Hzだと収束していないが，10Hzで収束
Structure A
Structure B

15Hz model
15Hz - 10Hz
15Hz - 5Hz

Structure B displacement EW (m)N UE• 形状を厳密にモデル化する手法を十分な分解能で使うことで，解が収束• 京＋GAMERAによって初めて都市ス
ケールで実現

15Hz model
15Hz - 10Hz
15Hz - 5Hz
27

ここまでのまとめ• 地盤震動問題は地震の統合シミュレーションの実現に際して，重要かつチャレンジングな問題– 構造物の被害を左右する問題だが計算量が膨大– Unstructured FEMの大規模・高速解析が必要に• 地盤震動問題のためのunstructured FEMの大規模解析手法を開発– 計算・通信を減らすことで，スケーラブル・高速な計算を実現– 京で実行することで，地盤震動の実問題を有用な分解能で解けるように– 形状を正確にモデル化する非構造格子を使うことで，綺麗に数理問題が解けている28

地震・津波問題のもう一つの側面• 実問題に適用する際には，地殻・地盤・構造物群などのモデルデータの欠如・不確実性も問題に– 例：広域の都市モデル．都市のデータは外形の3D情報はあっても，解析に使う構造データは整備されていない• 数理問題を高詳細に解くのと同時に，現実問題を計算に反映させる「モデル構築」も重要– データを有効活用する「自動モデル構築」と，– 不確実性を結果に反映するための「capacity computing」を合わせて結果の信頼性向上を試みる

29

Utilization of GIS data of buildings*
CAD data is just a set of polygons
independently located for visualization
Major Challenge:• Conversion of CAD data to models
with element connectivity informationObjective:• Develop a method for automaticconstruction of building models
Input: CAD data of buildings*Hideyuki O-tani et al.: Template-Based
Automatic Construction of Building
Models from GIS data, COMPSAFE 2014

CAD data (Polygons)

Slab-Column Model30

Automated Model Construction
horizontal polygons
with height

Floor arrangement

Floor with structure memberslab with beam

polygon 1
columnpolygon 2
wallpolygon 3

virtual polygon

CAD data

Floor arrangement setStructure Model31

Automated Model Construction• Two Major Challenges:• Floor shape recognition• Assumption of floor arrangementFloor shaperecognition

Assumption

columnwall

CAD data (polygons)

Floor shape

Floor arrangement32

Template Fitting
Floor arrangement

Floor polygons

convertTemplate framewith Floor arrangement33

Results of template fitting
(1,1)

(a,b)

(a,b)

(0,0)

rectangle(c,d)

L-shapeconvexThe results of template fitting applied to about 38,000 CAD data of buildings

Rectangle

L shapeConvex

Fail

0.05

47.8%

11.3%

1.3%

39.6%

0.10

54.6%

26.7%

5.9%

12.8%
34

Results of template fitting

d = 0.000

d = 0.05

d = 0.001

d = 0.05

d = 0.011

d = 0.335

d = 0.05

d = 0.565
35

Application

CAD Data (Polygons)

GIS Data (Polygons)

Slab-Column Model

Slab-Column Model

Slabやcolumnの位置を床形状に応じて推定可能36

構造物のstochastic modeling• 都市のモデル化に際しデータが欠如している場合がある– 例： コンクリートの物性．同じ設計でも，実際の性能はばらつく• このような不確実性も被害想定に反映させることが望ましい– Capacity computingによる多数ケースの解析で，入力分布に対する結果の分布を求める

Black dots: observed data by Prof. Noguchi Hirofumi,
Colored distribution: fitted bi-normal distribution

37

構造物のstochastic simulation• 13,275棟の都市を10,000サンプル計算– 京コンピュータ10,000ノード, 3h56minで計算– 応答の確率分布を求めることができるように
Structure A
Structure B

N UE38構造物のstochastic simulation• 全構造物に対して，応答分布を求めることができる– 対策の優先度を決める際などに有効と期待

層間変位角の平均応答・最大応答・標準偏差
39

おわりに
• 地震・津波問題は「広域・高分解能な物理」と，「不確実性を持つ大規模モデル」という二つの側面がある– 物理を高詳細で追うcapability computing– データを有効活用するためのautomated modelingと，その限界を補うcapacity computing (不確実性モデリング)– HPCはこれら２つの側面に対して重要な役割を担っている• ポスト京への期待– 地盤震動シミュレーションは京を使ってverifyできても，validationには足りない– ポスト京の計算能力を生かした多数シナリオの計算によるvalidationやstochastic simulationを行うことで，本当に役に立つツールとなることを期待40